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Abstract 

Using the well-known trigonometric six- vertex solution of the Yang-Baxter equa- 
tion we derive an integrable pairing Hamiltonian with anyonic degrees of freedom. 
The exact algebraic Bethe ansatz solution is obtained using standard techniques. 
From this model we obtain several limiting models, including the pairing Hamilto- 
nian with p + ip-w&ve symmetry. An in-depth study of the p + ip model is then 
undertaken, including a mean-field analysis, analytic and numerical solution of the 
Bethe ansatz equations, and an investigation of the topological properties of the 
ground-state wavefunction. Our main result is that the ground-state phase diagram 
of the p + ip model consists of three phases. There is the known boundary line with 
gapless excitations that occurs for vanishing chemical potential, separating the topo- 
logically trivial strong pairing phase and the topologically non-trivial weak pairing 
phase. We argue that a second boundary line exists separating the weak pairing 
phase from a topologically trivial weak coupling BCS phase, which includes the 
Fermi sea in the limit of zero coupling. The ground state on this second boundary 
line is the Moore-Read state. 
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1 Introduction 



Exactly solvable many-body models have played a very important role in various branches 
of physics, especially in statistical mechanics and condensed matter. They have been 
fundamental in the understanding of phase transitions and non-perturbative phenomena 
in low-dimensional and strongly correlated systems. Prominent examples are provided 
by the 2D Ising model, the ID anisotropic Heisenberg model, the interacting ID Bose 
gas, and the ID Hubbard model [H HI]. The solvability and integrability of these and 
related models emerges from a common algebraic structure, due to the quantum Yang- 
Baxter equation satisfied by the so called i?-matrix. The i?-matrix is the basic object 
of the Quantum Inverse Scattering Method (QISM) [11. It allows for the construction 
of a one parameter family of commuting transfer matrices whose expansion, in powers 
of the spectral parameter, generates an infinite set of conserved quantities including the 
Hamiltonian of the system. For the anisotropic Heisenberg model, also known as the XXZ 
model, the i?-matrix is the six-vertex trigonometric solution of the Yang-Baxter equation 
and it depends on a spectral parameter u and the quantum parameter q related to the 
anisotropy as A = (q 2 + q~ 2 )/2. 

A closely related family of exactly solvable models are the Richardson- Gaudin models, 
where the most prominent examples are the Richardson model of BCS superconductivity 
0, IIHI [Ml E2J EH] and the Gaudin spin models [23] (for a review see [2]). The common 
algebraic structure of this family is the classical Yang-Baxter equation satisfied by the 
classical r-matrix. This latter equation is obtained from the quasi-classical limit of the 
quantum Yang-Baxter equation through R(u) ~ 1 + hr(u). A feature of these models, 
unlike the anisotropic Heisenberg model, is the dependence of the Hamiltonian on a large 
number of parameters. In the Richardson model, they are the energies of the single 
particle levels and the BCS coupling constant, while in the Gaudin model they are the 
position of the spins. This abundance of parameters has its origin in the transfer matrix 
which is the product of i?-matrices with shifted spectral parameters. This is called an 
inhomogeneous transfer matrix. This situation must be compared with the anisotropic 
Heisenberg model whose transfer matrix is homogeneous, which is the reason why the 
Hamiltonian of this model only depends on q. 

Thus the majority of models constructed so far are based either on a homogeneous 
transfer matrix (e.g. anisotropic Heisenberg model), or in the quasi-classical limit of 
an inhomogeneous transfer matrix (Richardson-Gaudin models). It is also of interest 
to construct physical models from an inhomogenous transfer matrix before one takes 
the quasi-classical limit. The first result of this paper is to construct such a model 
starting from the trigonometric six- vertex i?-matrix discussed above. This model describes 
the pairing interaction of anyonic particles with braiding statistics in momentum space 
parameterized by g, and for that reason it can be termed an anyonic pairing model. These 
anyons should be distinguished from the real space anyons which satisfy braiding statistics 
in ID |] or in 2D HUES]. 

The one-body version of the anyonic pairing Hamiltonian we obtain coincides with 
the quantum mechanical Hamiltonian proposed by Glazek and Wilson in 2002, which 
furnishes a simple example of a renormalization group with limit cycles [2U [25] . The 
latter work inspired the construction of a BCS model with s-wave symmetry and time 
reversal symmetry-breaking where the renormalization group also has limit cycles [15] . 
The mean-field gap equation of this BCS model has an infinite number of solutions related 
by discrete scale invariance, which explains its name: Russian doll BCS model. This model 
was shown in reference [16] to be exactly solvable by the QISM using a rational R matrix. 
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Figure 1: Web of integrablc pairing models. 



The relation between the mean-field and the exact solution of this model was established 
in reference [I]. It then comes as no surprise that the Russian doll BCS model can be 
obtained from the anyonic pairing model in the limit where the trigonometric i?-matrix 
becomes rational. However, if the quasi-classical limit of the anyonic pairing model is 
taken before the rational limit a BCS model with p + zp-wave symmetry is obtained, as 
has been reported in [53] • Taking a rational limit of the latter model or a quasi-classical 
limit of the Russian doll model one arrives to the well known s-wave BCS model of 
Richardson [51]. Finally, one can take the limit where the BCS coupling constant, of both 
the s- wave and p + zp-wave models, goes to infinity obtaining two Gaudin models. The 
web of relations between all these models is shown in Fig. <^f. 

The first part of this paper, Sect. 2, is devoted to the construction of the anyonic 
pairing model using the QISM and the derivation of the eigenvalues of the transfer matrix. 
In Sect. 3 the relationships of the models in Fig. [T]are shown, giving a unified picture of 
the pairing models that can be obtained from the trigonometric R— matrix. In Appendix 
A we collect the more technical computation of the wavefunction and scalar products of 
the anyonic model, as well as its relation with the aforementioned Glazek- Wilson model. 

The bulk of the paper, Sect. 4, is devoted to a detailed study of the exactly solvable 
p + ip-w&ve model which was initiated in [33J . In recent years this model has received 
considerable attention due to the existence of vortices with non-trivial braiding statistics 
[53] |32| 160] . with the potential for implementing topological quantum computation [50] . 
Furthermore, there are some experimental indications that Sr 2 Ru04 : is such a chiral p- 
wave superconductor [12] , supporting topological vortices [121 [33] . In the realm of cold 
atoms there are also possibilities to realize p + zp-wave paired phases (TTJ [211 EU 03J 
through experiments involving Feshbach resonaces [211 1221 [3U EH EH [67] . 

We begin our study of the p + ip model using the BCS mean-field theory, to determine 
the phase diagram in the plane (l/g,x), where 1/g is the inverse of the BCS coupling 
constant and x is the fraction of occupied energy levels. This phase diagram is divided 
into three regions denoted as weak coupling, weak pairing and strong pairing in the 
terminology introduced by Read and Green [53]. The strong and weak pairing regions are 
separated by a line in the (1/g, x) plane where the chemical potential vanishes, giving rise 
to a second-order phase transition and a topological transition described by a winding 
number in momentum space introduced by Volovik [50] • The latter two regions are related 
by a strong-weak duality which connects two ground states having the same total energy, 
and value of the gap, while the chemical potential differs in the sign. On the other hand, 
the weak pairing and weak coupling regions are separated by another line in the (1/g, x) 
plane which can be called the Moore-Read line because the projected BCS state coincides 
precisely with the Moore- Read pfaffian state related to the Fractional Quantum Hall state 
at filling fration 5/2. All the points of this line have zero total energy and are mapped, 
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under the duality transformation, into the vaccum state. However the gap does not vanish 
on this line and the nature of the transition is rather subtle since, as we show below, it 
is associated to a discontinuity of the ground state energy, or in another words a zeroth 
order phase transition. Volovik has discussed the existence of a Higgs like phase transition 
for the p + ip model, meaning a qualitative change in the dispersion relation of the quasi- 
particles $T\. We find that this type of transition occurs on a line located in the weak 
coupling region, which is close, but not identical, to the Moore-Read line. A peculiarity 
of the Read-Green and the Moore-Read lines is that they are insensitive to the particular 
choice of the energy levels, and in that sense they are topological. This feature is not 
shared by the Volovik line whose position in the phase diagram depends on the energy 
level distribution. 

The analysis of Sect. 4 continues with the computation of the fidelity susceptibility 
which confirms, as expected, the singular nature of the Read-Green line, while the Moore- 
Read line does not show any special feature. Next we analyze the continuum limit of the 
Bethe ansatz equations of the p + ip model using an electrostatic analogy which allows 
us to uncover the complex structure of the arcs formed by the roots of the Bethe ansatz 
equations. We also recover the mean-field gap and chemical potential equations as was 
done in previous works for the BCS model with s-wave symmetry [231 ESI EH]. The 
structure of the arcs, the values of the ground state energy, gaps, chemical potential and 
occupation numbers are compared with the numerical results obtained by solving the 
Bethe ansatz equations for systems with finite size, finding a resonable agreement. A 
detailed analysis of the exact solution allow us to understand the weak-strong duality 
in terms of a dressing operation that maps exact eigenstates between those two regions. 
The physical picture that emerges is that the ground state of the weak pairing region is 
formed by strongly localized pairs and by delocalized pairs with zero energy. We then 
discuss the zeroth order phase transition and the winding number of the mean-field and 
exact solutions. For the mean-field solution the winding number vanishes in the strong 
pairing region while it does not in the weak pairing and weak coupling regions. However 
in the exact solution the winding number only vanishes in the weak pairing region where 
it counts the number of Bethe ansatz equation roots that are equal to zero. In this sense 
the winding number plays the role of a topological order parameter for the weak pairing 
region. Another qualitative distinction between the weak pairing and weak coupling 
regions is derived by studying the structure of the vortex wavefunction by means of the 
Bogolioubov-de Gennes equations. We find that the vortex wavefunction has an oscillatory 
behaviour in the weak coupling region which does not appear in the weak pairing region. 

Finally, in Appendices A and B we discuss several technical issues: relation between 
the Glazek- Wilson model and the anyonic model, computation of wavefunction scalar 
products and correlation functions of the anyonic pairing and the p + ip models, a q- 
deformed version of the Moore-Read state, solution of the gap and chemical potential 
equations. 

2 Construction of an integrable anyonic pairing model 

2.1 The Yang— Baxter equation 

A systematic method to construct integrable quantum systems is through the Quantum 
Inverse Scattering Method (QISM) [H]. The key element in this approach is a solution 
of the Yang-Baxter equation [H [38] , commonly known as an i?-matrix. This is a linear 
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operator acting on a two-fold tensor product space which is dependent oniGC known as 
the spectral parameter. Denoting the i?-matrix as R(x) G End(V (g) V), the Yang-Baxter 
equation reads 



R 12 (x/y)R 13 (x)R 23 (y) = R 23 (y)R 13 (x)R 12 (x/y) 



which acts on the three-fold space V £g> V <g> V. The subscripts above refer to the spaces 
on which the operators act. E.g. 

R 12 (x) = R(x) <g> /. 

The specific example we will use is the well-known XXZ solution associated with the 
quantum algebra U q 



( (q 2 x — q 2 x x ) 



R{x) 



(x — x x ) (g 2 



(f 








\ 



V 








(q 2 ~ q~ 2 ) 




(x — x x ) 




(q 2 x — q 2 x x ) J 

We define the local L-operator as L(x) = R{xq~ l ) which is a solution of 

R 12 (x/y)Li 3 (x)L 23 (y) = L 23 (y)Li 3 (x)R 12 (x/y) (2) 

as a consequence of ([T|). It is customary to consider L(x) G End(K ® V) as acting on the 
tensor product of an "auxiliary" space and a local "quantum" space. We use the subscript 
a to denote the auxiliary space and the labels j = 1, L for the local quantum spaces of 
a many-body system. 

To apply the QISM to construct a pairing Hamiltonian we next introduce the hard- 
core boson (Cooper pair) operators in terms of fermion operators. Let {Cj a , Cj a : j = 
1, L,a — ±}, denote a set of creation and annihilation operators with the usual canon- 
ical anticommutation relations and where the labels ± refer to time-reversed pairs. The 
hard-core boson operators are defined 



We also set Nj = bjbj as the Cooper pair number operator such that the total number of 
Cooper pairs in a system is given by the eigenvalue of 



N 



3=1 



Throughout we will use M to denote the eigenvalues of N. 

For each j the local Hilbert space of states is four-dimensional. Excluding the subspace 
of states which contain unpaired fermions reduces the local Hilbert spaces to dimension 
two, spanned by the vacuum |0) and the paired state bj\0). Now expressing the linear 
operators of the quantum space of the L-operator in terms of hard-core boson operators 
we have 



L a j {x) 



x 



q 



2N, 







q 



o 

I-2N, 



(q 2 - q- 



&„■ 



'J 







— x 



q 



I-2N, 







q 



o 

2Nj-I 



1 It is convenient for our purposes to express the deformation parameter as q 2 rather than the more 
familiar q. 
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noting that 

qt*i = (gT - l)Nj + I. 

Using the L-operator in this form, we will next apply the QISM to construct an integrable 
pairing Hamiltonian with anyonic degrees of freedom. 



2.2 Transfer matrix, Hamiltonian, and exact solution 

The most common use of the QISM is to construct a transfer matrix on a closed chain 
(to adopt the language of spin models). Here we use this approach with a "twist" in 
the boundary conditions, which will be chosen in a sector dependent manner as described 
below. The quantum space for such a model is the L-fold tensor product space W = V® L . 
We introduce the monodromy matrix T(x) G End(V ® W) defined by 

T(x) = U a L al (xZi 1 )L a2 (xz2 1 )....L aL (xzl 1 ) 

where 



U 



exp(—ia ) 
exp(ia / 

The monodromy matrix is commonly expressed in a matrix form as 



T(x) 



T±{x) T 2 V) 
T?(x) Ti(x) 



where the entries are operators acting on W. As a result of 02]) and 

[U ® U, R(x)} = 

it follows that 

R ab {x/y)T a (x)T b (y) =T b (y)T a (x)R ab {x/y) (3) 

which is an operator equation on V ® V <S> W, with the two auxiliary spaces labelled by 
a and b. 

Defining the transfer matrix to be 

t(x) = tr (T(x)) = Tl(x) + T|(x). 
it follows from ([3]) that the transfer matrices form a commutative family; i.e. 

[t(x), t{y)} = Vx, y e C. (4) 
The transfer matrices can be expanded in a Laurent series 



t (x) = t {j) 

3=~L 



X J 



and because of 01]) the co-efficients commute 

t {k) ] =0, —L < j, k<L. 
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One may check directly that the transfer matrix commutes with N, so the Hilbert space 
of states can be decomposed into sectors labelled by the eigenvalues M of N. The leading 
order terms in the expansion are 



t (L) 
t (L-l) 



exp(2«')g i " 27V ) 



0. 



We will construct a Hamiltonian in terms of the operator t^ L 2 ** , where the choice of the 
variable a will depend on the sector labelled by M. Setting 



q 

a 



exp(i/3) 



a 



we define 



H 



4sin(2/3) sin(a) 



f3{L-2M + 2) 



^ 2 )n^ + 2cos(«)^ 



j"=i 



The Hamiltonian may be expressed as 

L 



H 



sin(2/3) 
sin(a — 2/3) 



^ z k z r (exp(—ia)dld k + h.c.) 



with 



(5) 



di 



h n * mk - 2i > 
k=j+i 



(6) 



a, /3 are assumed to be real and h.c. denotes hermitian conjugate. Due to the non-local 
action of the operators dj, dj they satisfy 



djdj, 
d\d k 



q A d k dj 



'1 



l d k d\ 



j>k 
j>k 



amongst other relations, and therefore can be viewed as anyonic creation and annihilation 
operators. By the nature of the above construction, the anyonic Hamiltonian is integrable 
as it commutes with all terms in the expansion of the transfer matrix. That is 



[H, t®] = 0, 



J 



—L, 



establishing that the provide a set of conserved operators^. 

Recall that the derivation of the Hamiltonian was made using a restricted subspace of 
the Hilbert space obtained by excluding unpaired states. By noting that the second term 



2 In principle for integr ability to hold one would need to check that the number of independent con- 
served operators is not less that the number of degrees of freedom. 
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in (j5J), which describes the scattering of Cooper pairs, vanishes on unpaired stated it is 
straightforward to extend the action of the Hamiltonian to the full 4 L -dimensional Hilbert 
space. We define anyonic creation and annihilation operators {cij±, CLj± '■ j = 1,...,L} 
which are defined in terms of the canonical fermion operators through 



k=j+l 



n <r 

k=j+i 

where rij a = c^ a Cj a . These operators satisfy the relations 

{a j(T , a jp } = {a j+ , a k -} = 0, (7) 

{a ja , a) p } = 5 ap I, (8) 

a ja a ka = -qa ka a ja j > k (9) 

a] a a k a = -q^a^a^ j > k (10) 

and those relations obtained by taking hermitian conjugates. The usual fermionic com- 
mutation relations are recovered in the limit q — > 1. Then the Hamiltonian 

1 v 1 ^ 2 i \ sin(2/3) ^ / . + + \ 

H = -^2j [n j+ + nj_) - — - — z k z T [exp(-ia)a , r+ a , r _a k _a k+ + h.c.J 
j=i sm{a p) fe> ^ 

is an extension to the 4 L -dimensional Hilbert space which has the same action as ([5]) when 
the unpaired states are excluded. 

Having constructed the Hamiltonian through the QISM it is a standard calculation to 
obtain the exact solution through the algebraic Bethe ansatz. We simply present the key 
results. Explicit details relevant to the derivation of the formulae below may be found in 

[2H1 SQl SH E2]. 

The first step is to identify that the vacuum state |0) admits the properties 

Tl(x)\0) = a(x)\0), 

T}{x)\0) = 0, 

7?(x)|0> ^ 0, 

T 2 2 (z)|0> = d(x)\0) 



where 



a(x) = exp(-ia') l xz k 1 - qx 1 z k ), 



k=i 

L 



d(x) = exp(ia') J^J(gx2 fc 1 — g 1 x 1 



z k ). 



k=l 



We then look for eigenstates of t(x) in the form 

M 

iw>=n i ?(vi) i°) 



3 In standard BCS theory this phenomenon is well-known and goes by the name of the blocking effect, 
e.g. see [63]). 
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where Y = {yj}- Note that 

[T?{y j ),T?(Vk)]=0 

so the operators may be ordered arbitrarily. Using the algebraic relations amongst the 
Tj(x) operators, this leads to the eigenvalues A(x) of the transfer matrix being given by 

m\ goi sb m 



L -^22—22 

A(x) = exp(— ia ) I I (q~ xz^ — qx~ Zf.) I I 
fe=i 



i= i ^ " ^ 



L 



Af — 2^,2 2 2 

+ exp(za') ^(gxz^ 1 - q^x^Zk) 
fc=i 



x 2 — y 2 



exp(— ia 



fc=i 



I 



M 2 -4 2 

x — q Uj 



3=1 



a; 2 — y 2 



L M 2 4 2 

— q y 

+ exp(za) ^(x^ 1 - g~ 2 x _1 z fc ) JJ ■ 
fc=i j=i 



x 2 — y 2 



such that the parameters {yj} satisfy the Bethe ansatz equations 



exp(-2za) JJ - 



L 9 9 9 M 

y m 4 

- 2iji — 2 ^2 
fc=l " " J ym fc jym 



q~ 4 y™ 2 yj 



m = 1, ...,M. 



'ill 



To obtain the energy expression i? for the Hamiltonian 
the transfer matrix eigenvalues in the Laurent series 



it is a matter of expanding 



giving 



A(z) = A ° )xj 

j=-L 



E 



4sin(2/3) sin(ot) 



3=1 3=1 



sin a 



Al 



sin 



:« - 2/3) 



(12) 



We leave to Appendix IA.2I the computation of the scalar product and correlators for this 
model. 

From the form of the Hamiltonian ([S]) and the definition of the operators dj, eq. (jH]) 
one can see that the model parameters can be restricted to the domain < a, (3 < it. If 
a = the Hamiltonian (J5J) reduces to 



H = z T d\ ^ z kdk 



r=l 



,k=l 



Formally eq. f ll2p yields that the energies are zero. However in this limit the roots yj can 
be divergent, such that non-zero energies do occur. 
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It is possible to generalise the above construction in a number of ways. Instead of 
using the closed chain transfer matrix construction, one can adopt the open chain version 
developed by Sklyanin [58]. In this manner a Hamiltonian is obtained which not only 
involves scattering interactions between the anyonic pair operators given by ([6]), but 
also their complex conjugates. A second route is to use different solutions of the Yang- 
Baxter equation. An obvious choice would be to employ Baxter's eight-vertex solution 
[7J |S]. But this solution suffers from the fact that it is not U(l) invariant which prohibits 
the construction of a Hamiltonian which conserves particle number. Beyond this there 
are however many known solutions of the Yang-Baxter equation which do possess U(l) 
symmetries, particularly those associated with representations of the quantum algebras 
U q [g\ where g is a classical simple Lie algebra (e.g. see [37]). In principle these can be 
applied for the construction of models, which in particular cases will provide anyonic 
generalisations of integrable pairing models based on bosonic or fermionic degrees of 
freedom. 



3 Limiting cases 

3.1 Russian doll and s-wave BCS models 

An example of a many-body system which admits a cyclic Renormalization Group (RG) 
equation analogous to the Glazek- Wilson Hamiltonian (see f 1 1 3 j) of Appendix lA.lj) was 
studied in [15]. The many-body Hamiltonian was taken to be the s-wave BCS model with 
complex-valued coupling parameter, later to become known as the Russian Doll (RD) 
BCS model. The cyclic RG was discovered through a mean-field analysis. Subsequently 
in [16] it was proved that the RD BCS model is integrable, and in [1] the connection 
between the Bethe ansatz solution and the RG equation was exposed. Next we will show 
that the RD BCS model can be obtained as a limiting case of the anyonic pairing model 
05]). Moreover the exact solution for the usual s-wave BCS model, which was first found 
in [51], is obtained by performing a second limiting procedure. 

We introduce a new variable 77 and define the parameters Ef., k = 1, L and Vj, j = 
1, M through 

z k = exp(2(3e k /7]) 
Vj = exp(2 (3vj/r]) 

and redefine the Hamiltonian (J5J) by a simple rescaling and additional conserved term 

— ^—(H-N) 
sin(2/3) v J 

so the energy expression flT2|) becomes 

• f \ M 
_ 7? sin (cm 9 V 

E = — > v ' — M. 

sin(a - 2/3) sin (2/3) yj sin(2/3) 

Now we take what is known as the rational limit (3 — > to obtain the following expressions 
for the Hamiltonian, Bethe ansatz equations, and energy respectively 

L L 

H = 2^e j N j -G^{e^{-ia)blb k + h.c) (13) 

j=l k>r 
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exp(-2za)|l — — = I — , j = l,...,N 14 

Af 

£ = 2^^- + GMcos(a). (15) 

3=1 

where G = r]/ sin(a). This is the exact solution for the Russian Doll BCS model in 
terms of the parameterisation used in [16J. Applying the same limiting procedure to the 
expressions ( I112|117p derived in Appendix IA. 21 yields the wavefunction scalar product and 
one-point functions given in [To] . 

Performing a second limiting procedure to the Russian Doll Hamiltonian yields the 
s-wave Hamiltonian. To make this result explicit we implement the following change of 
variable a — > rja and add the conserved quantity —GN to the Hamiltonian ( 1T5|) . Taking 
the quasi- classical limit t] — >• we obtain from (I13|ll4|ll5p 

L L 

H = 2^e j N j -G^b% (16) 

j=l k,r=l 



M 

V 



= > j = h~,M (17) 



7 1 3 ^ k --L 

fc=l J J7^ m 

M 

E = 2^2 



V 3 



3=1 



where G = 1/a. Up to a change in notational conventions this is precisely the exact 
solution first given in [54] . 



3.2 The p + ip-wave BCS model 

The Hamiltonian of the BCS pairing model with p + zp-wave symmetry was introduced 
in reference [33] 

k 2 G 

H = E 2^ C ^ Ck ~ 4^ E ~ *^)(^ + ik 'y) c k c -k c -k'C k ' (19) 

k k^±k' 

where destruction and creation operators of two-dimensional polarised fermions 

with momentum k = (k x , k y ), m is their mass and G is a dimensionless coupling constant 
which is positive for an attractive interaction. We remark that we impose no constraint 
on the choice for the ultraviolet cut-off, which we denote as u. Likewise the distribution 
of the momenta k is arbitrary other than the assumption that all momentum states arise 
in time-reversed pairs k and — k. In particular this means that a one-dimensional system 
is obtained by simply setting all k y = 0. 

We shall discuss now how this Hamiltonian can be obtained from another limiting 
case of the Hamiltonian (j5J), which establishes the integrablity of the p + ip-wave BCS 
model. The strategy here is to take the quasi-classical limit first, rather than the rational 
limit. While this general approach has previously appeared in the literature to construct 
integrable Hamiltonians in the context of Gaudin algebras [U [2j [TH [15] , the connection 
with the p + ip model has thus far not been made apparent for reasons we will discuss 
below. 
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Introducing the parameterisation 

a = —ijt 
(3 = —i'jp 

and taking the limit 7 — > leads to the Hamiltonian 



j=l k>r 



with the Bethe ansatz equations and energies given respectively by 

l 2 m 2 

-^-^E-T^ = -^EA' ( 2 °) 

k=l Um k j^m y j y ™ 

t M 

E = T^X* 



Next we set 



in terms of which we have 



t-2p a/fi-2 ( J 



j=l k>r 



G~ 1 -L + 2M-1 A 1 A 2 



M 

E = (i + g)£# ( 24 ) 

To show that (|22|) is equivalent to (fl9|) . we first define the Cooper pair operators 

= c k c -k> &k = c_ k c k , iV k = blb k , 

which have odd symmetry (consistent with p-wave pairing): 

6-k = -6k, iv k = AL k . 

We let K + denote the set of momenta where k x > and no restriction placed on k y , so 
that we avoid issues with double counting. Excluding the unpaired states and setting 
m — 1, the Hamiltonian f fT9|) takes the form 

H = 22 zl Ni, - G 22 ZkZk ' ex P(-^0k)exp(i0 k /)6j c 6 k / (25) 

k6K + k^k'GK+ 

where z k = |k| and exp(i0 k ) = (k x + ik y )/h 2 . Performing the unitary transformation 

b\. = exp(i0 k )6j c , b k = exp(-i0 k )6 k (26) 
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and using the integers rather than k to enumerate the momentum states in the half plane 
brings Q25J) into ffT9|) . 

Note that the main results derived for the exact solution of the p + ip model, viz. 
(I23|24p . are functions of the squares of the Bethe roots. Throughout the subsequent 
analyses (with the exception of Appendices IB. 31 and IB.4p we will simplify the notation by 
making the substitution 

V) ^ Vr (27) 

In terms of the original parameterisation, we may now write that the exact eigenstates 
of the Hamiltonian f lT^]) with M fermion pairs are given by 

m = 11^)10), C(y) = C U (28) 

J'=l k6K+ * 

where the rapidities yj, j = 1, M satisfy the Bethe ansatz equations 
1 i M i 

with 2g= 1/G — L + 2M — 1. The total energy of the state (128]) is given by 

M 

J5 = (l + G)5>. (30) 

This limiting case is closely related to other models which have been constructed using 
the hyperbolic Gaudin algebras [H [21 [HI [15] . In fact the conserved operators for the 
Hamiltonian f[22|) which are obtained by taking the quasi-classical limit of the operators 
{t(qzj) : j = 1, ...,L} are equivalent to those which are discussed in [H [21 [TU [15] with 
an important difference due to the coupling parameter of the overarching anyonic model 
being defined in a sector-dependent manner. To make this explicit, expressing the leading 
term expansion as 

t(qz k ) ~ / - 2p'jr k 

we find 

r k = (2G- 1 -L + 2M)N k - J2 (4^2 (t>lh + b k b\) + 4^4(2iV fc iV, - N k - .(31) 

;_^ fc \ Z k Z l V 7 % Z « / 

Besides a constant term, the difference with the conserved operators in (TJ [2} H3J [15] is 
the M- and L-dependence on the co-efficient of iV^. Through the conserved operators we 
can write for each sector with fixed M 

G 



k=l 

The eigenvalues of the r k are 



M 2,2 
\ ^ Z k + ^ 
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which follows from taking the leading terms in the expansion of f ll09p in Appendix IA.2I 
One may check, with the help of (123"]) . that 



E = (1 + G)£>? 



M 



L 



fe=l 

as required. 

It is important to clarify here the reason why the conserved operators are said to be 
associated with hyperbolic Gaudin algebras. Introducing the scaling parameter v we set 

e fc = -ln(zfc), Vj = -\n(yj). (32) 

Through this change of variables the Bethe ansatz equations (1231) can be expressed as 

L M 

_ t, 4- 9.M + " 

G 



^ -L + 2M + s ^coth{v{v m -e k )) = 2 E coth(z/(t; m - Vj )) (33) 

fc=l 

and in a similar way the expressions for r k and X k can be written in terms of hyperbolic 
functions. In this form, we can now take the rational limit to recover the s-wave BCS 
model from the p + ip model. We redefine the Hamiltonian ( 1221) as 

H h-> -(H- (1 + G)N). 
v 

Then setting G i— > vG and taking the limit v — > 0, ( )22|23|24p reproduce ( JT6|17|18p . 



3.3 The Gaudin models 

The Gaudin model can be considered as the G — > oo limit of the corresponding s- wave and 
p-w&ve models introduced above. In the s-wave model the U(l) symmetry is promoted to 
a SU (2) symmetry and the Gaudin model can be regarded as a rotational invariant spin 
chain with local spin 1/2 operators Sj whose dynamics is described by a set of commuting 
operators 

Ri = ^ j , [Ri, Rj] = 

These operators are simultaneously diagonalized by the Richardson ansatz in terms of a 
set of rapidities e,j satisfying the G — > 00 limit of the Bethe ansatz equations (fTTl) . namely 

L M 

^— ' t> j — e fc f m — v.- 

k=l J j^pm 

The p-wave Gaudin model corresponds to the trigonometric version considered earlier 
with the commuting operators given by the G — > 00 limit of the operators r k defined in 
eq. (I3"T|) . The corresponding Bethe ansatz equations are given by the G — > 00 limit of 
or equivalently (1531 . 
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4 Analysis of the p + ip model 



4.1 Ground-state phase diagram from BCS mean-field theory 

To undertake the mean-field analysis of the p + ip model it is more convenient to work 
with the Hamiltonian 

k 2 G 

U = Yl 2^ C k Ck ~ 4^ J2( kx ~ ik v)(K + ik 'y) c k c -k c -k'C k ' (34) 

k k,k' 

which differs from f lT9|) only in that the double sum is now not restricted to k ^ ±k'. The 
differences between using (fT§j) and (154]) will be discussed at the end of Subsection 14.1.11 
The derivation of the gap and chemical potential equations which give the mean-field 
solution of (134]) follows now standard techniques which are a straightforward extension of 
the original methods of the BCS paper [5]. Here we will only outline the main steps. First 
we assume that the number of fermions is even and that in the ground state of the system 
all fermions are paired. Our convention will again be to use N to denote the Cooper pair 
number operator and M for the eigenvalues of N. The BCS order parameter associated 
to flS} is 

A = — y~)(k x + ik y )(c^ k c k ) 

k 

in terms of which the Hamiltonian (1541) can be approximated as 

^ - 7 ( A " iA v)i c -k + h - c ) + ^ + / iM 
k k 

where £k = k 2 /2m — /i/2 and fi/2 is the chemical potential. This Hamiltonian can be 
diagonalized by a Bogoliubov transformation. By minimising the energy, and fixing the 
expectation value of the number of Cooper pairs to be M, it is found that the gap A = |A| 
and chemical potential are the solutions of the equations 

k 2 1 

7 = 77 35 

ihU V(V - /i) 2 + k 2 A 2 G 

H V 1 = 2M-L + -7 (36) 

ktk^ + V(k 2 - /^) 2 + k 2 A 2 G 

where we have set m = 1 and 2L is the total number of momentum states below the 
cut-off w. The mean-field expression for the ground-state energy is 

The normalised ground state is given by 

l*> = II (^I + vAd k )\0) (38) 

keK+ 
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where 



11 k 2 - a 



2 V v/(k 2 - /i) 2 + k 2 A 



2 



k 2 -/i 

= " I 1 



2 \ v/(k 2 - /i) 2 + k 2 A 2 
and the phases must be chosen such that 

is real. Above, E(k) is the quasiparticle energy spectrum 

E(k) = i v / (k 2 -/i) 2 + k 2 A 2 . (39) 

To obtain an approximation for the ground state with a fixed number of Cooper pairs 
M, one can take a projection of fl38|) . Writing (155]) as 

i*> = ( n ^ ) ( e ) 

\keK + / \k 6K+ nk / 
projection onto a fixed number of M pairs gives 

-l M 



keKj 



|0> (40) 



where 



(k) = * = 2E(k)-k^ + ^ 

(k x + ik y )A* 

4.1.1 The Moore- Read and Read-Green lines 

By analysing the above mean-field results we can piece together the various phases of 
the model. To some extent this has been undertaken in [53], and we will first review 
their findings. Then we will extend the analysis to expose a duality that exists in the 
ground-state phase diagram. 

Following [53], note that the spectrum is gapless at \x = as |k| — > 0. Furthermore, 
the behaviour of g(k) as |k| — > depends on the sign of /i, 



0(k) 



k x iky, /i < 0, 
l/(k x + ik y ), /i > 0. 



It was argued in [53] that /i = corresponds to a topological phase transition, a subject 
we will return to later in Subsection 14.31 From ( 13 6 p it is seen that \x = necessarily 
implies 



2 \ G 
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Introducing the filling fraction x = M/L and defining g = GL we equivalently have 



2 V 9 

We will refer to ( 142]) as the Read-Green (RG) line of the phase diagram. 
In real space the state (140]) takes the form of a pfaffian 

ip(r ll . . . , r 2M ) = AQ( r i ~ r 2) • • • 0O2A/-1 - r 2M )] 

where A denotes the antisymmetrization of the positions and $j(r) is the Fourier transform 
of g(k). We will refer to the case /x = as the Read-Green (RG) state. For \i > the large 
distance behaviour is $j(r) ~ l/(x + iy), which asymptotically reproduces the Moore- Read 
(MR) state found in studies of the Fractional Quantum Hall Effect [49] (see Appendix 
IB. II for the derivation of this result and its comparison with the exact wavefunctions in 
real space). In momentum space co-ordinates the (unnormalised) MR state of M Cooper 
pairs is 



I Mi?) 



v — 

t-~< k... 4- 



-1 AI 



■ k C k C -k 



|0>. (43) 



Upon closer inspection of the mean-field equations we find that the MR state (1431) 
coincides (up to normalisation) with the projected ground state (HOI) whenever 

A 2 = 4/x (44) 

The result is easily verified by simply substituting f!44p into (I39II41I) . Furthermore, sub- 
stituting (jHJ) into (I35|I36I) and adding these equations yields 



or in terms of the filling fraction 



1 , X 

x = 1 - -. 45 

9 



We will refer to (145*]) as the MR line of the phase diagram. Finally, we substitute 
into f[37]) to obtain 

E = 0. (46) 

This means that the condensation energy of the state equals the energy of the Fermi sea! 

From the above considerations we gain an initial insight into topological properties of 
the model. Both the RG line given by (142 p . for which the excitation spectrum becomes 
gapless as |k| — > 0, and the MR line given by (1451) . for which the ground-state energy 
is zero, hold independent of the choice of distribution for the momenta k. These are 
properties of the model that are protected (in a topological sense to be discussed in 
Subsection 14. 3p from perturbations of the system which are reflected by changes in the 
momentum distribution. 

Both the RG and MR lines can be viewed as boundary lines in the ground-state phase 
diagram relative to a duality property of the mean-field solution, which is our next item 
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to discuss. First we introduce parameters a and b related to the chemical potential and 
the gap by 



A 2 /A 2 

a,b = fi- — ±A^—-fi (47) 



which simply allows us to write the quadratic term under the square root in ( )35|36l) in 

terms of its roots: 

V(k 2 - /if + k 2 A = v/(k 2 -a)(k 2 -6). 

The parameters a and b can either both be real, or form a complex pair in which case we 
denote them as e ± i8. In either case they satisfy the relation 

/j 2 = ab. 

We see from (157)1 that in terms of a and b the ground-state energy reads 

1 , , / 2k 2 -a-b \ 

E = - V k 2 1 == 48 

2^ + ^ 2 v /(k 2 -a)(k 2 -6 ) ; 

which is independent of the sign of \i. On the other hand, the sign of the fi coincides, by 
eq. (|36|) . with that of 

(49) 

This implies the existence of a duality between phases which we will denote as weak 
pairing and strong pairing, adopting the terminology of [53] . Specifically, the two ground 
states with filling fractions xw (with ji > 0) and x$ (with \i < 0) satisfying 

xw + xs = 1 - - (50) 

share the same values of a,b,g,fi 2 . Besides differing in the sign of /i, they also differ in 
the values of A. One can check that the required map which preserves (T47j) is 

A 2 ^ A 2 -A/2 

\i 1—7- —\l 

Finally, from (I48p we see these ground states have the same mean-field energy, and from 
03n]) the same excitation spectrum above them. 

Taking into account all these considerations, the ground state of the model consists of 
three phases detailed in Table 1. In the weak coupling BCS phase where x > 1 — 1/g, a 
and b form a complex pair parameterized by e and 5. When x = 1 — 1/g, the MR line, 
both a and b are equal and real. This line is dual to the vacuum. From ( 13 9 p the mean-field 
solution predicts this line is gapped, so it is not clear at this stage in what sense the MR 
line might be considered a phase transition. We shall come back to this issue later on 
when we discuss the solution of the model in terms of the Bethe ansatz equations. The 
two regions denoted weak pairing and strong pairing, for which the parameters a and b 
are real and negative, are separated by the critical RG line where the chemical potential 
/i — and the excitation gap vanishes. The RG line is self-dual. 
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Phase 


a, o 




X 


Weak coupling BLb 


e±iS 


\ '} I A 

fx > 


x > 1 — g 


Moore-Read line 


a = b = —fj, 


A '2 / A 

p = A z /4 


i -1 
ZJVffl = 1-5 


Weak pairing 


a < b < 


< /x < A 2 /4 


(l- ff - 1 )/2<x< l-^ 1 


Read- Green line 


a < b = 


/i = 


Xi ?G = (l-5" i )/2 


Strong pairing 


a < 6 < 


< 


x<(l- 5 - i )/2 



Table 1: Phases of the p + ip model in terms of the mean-held variables fi, A, a, b. 

To solve the mean-field gap ( 135]) and chemical potential fl36|) equations for the case 
of large L it is useful to undertake a continuum limit whereby the sums are replaced by 
integrals. Making the change of variable e = k. 2 /u we introduce a density function p(e) 
with the choice of normalisation 



de p(e) 



1. 



The densities are chosen as 



(51) 



which corresponds to free fermions in one dimension (ID), and 

p{e) = 1 (52) 
for free fermions in two dimensions (2D). The integral approximations of (I35|36p are 

ep(e) 



1 

9 





1 


1 




2 






= \f~ab 




X "2 H 


~% 


Jo 



de- 

o \/ (e — a){e — b) 
Pie) 



de- 



(53) 
(54) 



where a = a/oo, b = b/u. From ( 1371) we obtain the ground-state energy per pair as 

2e — a — b 



e 



2x 



deep(e) 1 



2^(e-d)(e-b) 



(55) 



The integral equation approximations become exact in the thermodynamic limit L — > oo, 
which requires that G — > 0, M — > oo in order for g = GL, x = M/L to remain finite. 
Solutions of ( I53|l5"4l) are given in Appendix IB. II Numerical results are presented in Figs. 
E] and H] below. 

To conclude here, we turn to the comparison between the Hamiltonians f fT9|) and fl34l) . 
Excluding the subspace of unpaired states we have 

H(G) = (1 + G)H{G/{1 + G)) 

which shows that H is equivalent to Ji up to a scaling factor 1 + G and the redefintion 
of the coupling G i— > G/(l + G). For systems with an even number of fermions, the 
mean-field ground states of H and % are equal in the thermodynamic limit since we take 
G^O. 
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4.1.2 The Volovik line 



The energy gap between the ground and first excited states is the minimum of the quasi- 
particle excitation energies as given by (139]) . In order to make comparison with later 
results for finite-sized systems where the total fermion number is fixed, we will consider 
an elementary excitation to be associated with two quasiparticle excitations. From ( 13 9 j) 
we find that the gap E gap is 

r h /i<A 2 /2, 

Egap(k) = { ^>/A2(4/i-A2), n > A 2 /2. 

For /i > A 2 /2 the quasiparticle energy has a minimum at a non zero momentum k 
different from zero, while for \i < A 2 /2 the minimum is achieved at k = 0. This situation 
is reminiscent of the Higgs transition that takes places for a scalar field with a quartic 
potential [61]. For this reason we call the case /i = A 2 /2 the Volovik line in the p+ip phase 
diagram. In this instance eliminating A from the gap and chemical potential equations 
( 1351136]) gives 



keK + 



1 



^ /k 4 + fi 2 G 



kGK+ 

Further eliminating the 1/G term we obtain the equation of the Volovik line 



which lies in the weak coupling BCS phase. In contrast to the RG (142]) and MR ( 145]) 
line equations, the Volovik line is dependent on the distribution of the momenta k and 
consequently is not topologically protected. In the continuum limit the equation ( 1561) 
reads 

x = -- - f 1 de p(e)^ZJL= (57) 

which relates the filling fraction x with the value of the normalized chemical potential fi. 
In 2D the gap equation (1531 becomes 

1 

9 



which can be easily inverted to yield /i as a function of g. Plugging that result into the 
integral of ( 1571) gives the Volovik line for the 2D model 



x = U\ - g- 1 + kg - g- 1 ) log 9 + 9 ~ 1 + 2 ) . (58 ) 

2 V. 2 g- g 1 J 



In ID the gap equation (1531) yields 

g 3/2 ' V 2 ' 4' 4' /i 2 
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Figure 2: Ground-state phase diagram of the p + ip model in terms of the inverse coupling 1/g and 
filling fraction x = M/L. Phase boundaries are given by the Read-Green line (/i = 0) and the Moore- 
Read line (fx = A 2 /4). The phase boundary conditions are independent of the choice of the momentum 
distribution, and independent of the ultraviolet cut-off. Also shown are the Volovik lines corresponding 
to ID and 2D momentum distributions. 

which together with the integration of fl57|) gives 

which is the equation of the Volovik line in ID. In Fig. [2] we depict the Volovik lines, of 
the ID and 2D models, on the the phase diagram of the p + ip model. One can see that 
they lie close to the MR line although do not coincide with it. As explained earlier the 
Volovik line depends on the choice of the energy level density p(s), while the MR and RG 
lines are independent of it. 

4.1.3 Fidelity susceptibility 

The notion of fidelity has been adopted in a variety of ways to characterise quantum 
phase transitions [651 EE]- Here we focus on one such formulation, which is the fidelity 
susceptibility [28] . 

The fidelity susceptibility \f is defined as 

dm \ 
dG / ' 

For the mean-field wavefunction this reduces to 




The derivatives of and with respect to G can be computed in terms of 

d/u dA 
~dG =C ' dG = ' 



Xf 



d^ 

dG 



d^ 

dG 



dm 

dG 



m ) ( m 



22 




1.4 
1.2 
1 

0.8 
fl 0.6 
0.4 
0.2 


-0.2 



2=1/4 — 
2=1/2 — 
- x=3/4 — 













0.5 



L.'i 



2 

.9 



(a) 



(b) 



Figure 3: Mean-field order parameter, A, and chemical potential, /i, solution of the gap and chemical 
potential equations in 2D for filling fractions x = 1/4, 1/2 and 3/4, vs. 5. The rectangles over the curves 
indicate, from left to right, the Volovik, Moore-Read and Read-Green points (x = 1/4); the Volovik and 
Moore-Read points (x = 1/2), and the Volovik point (x = 3/4). 




Figure 4: The mean-field energy gap E g&p versus coupling constant g for a 2D momentum density 
distribution as given by (|52j) . From top to bottom the curves shown correspond to the filling fractions 
x = 3/4, 1/2, 1/4. The squares over the curves indicate, from left to right, the Volovik, Moore-Read, and 
Read-Green points respectively. 



which never vanish for G > 0, as can be seen from Fig. [3j Now 



dui 



dG 



e k (cA + d(e k - n))' 



dG J ' \dG J 4 ((e k - /^y + e k A 2 ) 2 

In the mean-field theory gapless excitations occur when eo = fi = 0. Now 



lim lim 

e — >0 p.— >0 



duo 
~dG 



dv 
dG 



00. 



This indicates the fidelity susceptibility is divergent on the Read- Green line. The expres- 
sion of xf in the continuum limit is given by 
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Figure 5: Plot of the mean- field fidelity susceptibility (|5§)) as a function of the coupling g for x = 1/4. 
The 2D momentum distribution (|52j) is used and the fidelity susceptibility is in units of L 3 . The singularity 
at g c = 2 is apparent. Arrows indicate from left to right the Volovik, Moore-Read, and Read-Green points 
respectiveley. The inset emphasises the logarithmic character of this singularity, with the distance to the 
critical point g c shown on a logarithmic scale. The two data sets correspond to both cases g > g c and 
.9 < 9c- 



XF 



de p(e) e 



cA + d(e — /i) 
{e-p) 2 + eA 2 



(59) 



— = d. 



where 

dp, _ dA 
dg dg 

In Fig. Owe plot xf as a function of g for x = 1/4. The results show a logarithmic 
divergence at the Read-Green point g c = 2, reflecting that the model is critical here. 



4.2 Analysis of the Bethe ansatz equations in the continuum 
limit 

4.2.1 The electrostatic analogy 

For the s-wave BCS model, the relationship between the mean-field gap and chemical 
potential equations and the exact Bethe ansatz solution was first explored by Gaudin 
|23j . This was achieved by use of an electrostatic analogy, whereby the Bethe ansatz 
equations represent the equilibrium conditions for a two-dimensional system of fixed and 
mobile charges. Further studies along this line for the s-wave case can be found in [551 EH] . 
Our aim below is to extend this approach for the p + ip model. In the general framework 
of exactly solvable models based on hyperbolic Gaudin algebras this type of study has 
previously been conducted in [TJ. There the mapping from the hyperbolic form of the 
Bethe ansatz equation such as (|33|) to a set of algebraic equations was implemented by 
a change of variable. However this change of variable is not the same as ( 132]) . so the 
algebraic form considered in [TJ is different from ( I23|) . 

Taking into account the change of notational conventions as described by ( 12 7p . let us 
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write the Bethe ansatz equation f[2"31 as 



L M 

+ =«. m = l,... ) M (60) 

where go is given by 049|) . This equation admits an electrostatic analogy associating —1/2 
charges at the fixed positions z%, a —go charge at the origin, and mobile +1 charges with 
positions yj. With these assignments (160]) amounts to the vanishing of the total electric 
field acting on the charge y m . The continuum limit of (160!) is given by (we shall follow 
the conventions of references l23l [551): 



/• ie *L_*_ P / W |lM =0 , vyzr (6i) 

Jn e-y y J r y'-y 

where Q = (0,u) denotes the interval on the real line where lie the energy levels which 
can be associated to a negative charge density —p(s) such that 

dep{e) = -. 

The function r(y) denotes the charge density associated to the roots y m which lie on an 
arc T of the complex plane. It satisfies 



r 



\dy\ r(y) = M 

while the energy (|24p in the continuum limit is given by 

\dy\ y r(y) = E. (62) 



r 

For the ground-state roots the topology of the arc V depends crucially on the domain of 
the phase diagram one is exploring. In each phase V consists of the union of several types 
of arcs whose properties are given in Table 2. Table 3 displays the arcs associated to the 
different phases of the model. We shall first comment on the results contained in these 
tables and will give the derivation later. 

An arc of type Ta is a segment of the real line (0,£a) which is contained in the real 
interval Q. At G = 0, all the roots y m associated to the ground state lie on Ta, where 
Ea = p only depends on the filling fraction x. When the coupling is switched on, g > 0, 
the roots closer to the Fermi level form a complex arc Tb parameterized by its end points 
a, b = e ± i5, and which cuts the real axis at Ea whose value has also changed slightly. In 
the weak coupling BCS phase, T is the union of the arcs Ta and T B . This type of arc is 
the familiar one encoutered in the solution of the s-wave model [23l [55] . As g increases, 
the arc Tb enlarges and bends towards the negative real axis. At the Moore- Read line 
the complex arc closes and is denoted by Tc- This arc cuts the real axis at Ea > and 
a = b < 0. Some real roots may still remain, so that T at the Moore- Read line is the union 
Ta^Tq- In the weak pairing phase the closed arc Tq remains but there appear some new 
roots on the negative real axis forming the open arc To = (a,b). Hence T = Ta UTc UT^. 
As g increases further, the closed arc Tc shrinks and Tp becomes larger until g reaches 
the Read-Green line, where Tc disappears. Finally in the strong pairing phase there are 
only negative real roots. In Fig. O we display the types of arcs depending on the region 
of the phase diagram together with the exact numerical roots. 
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Figure 6: Representation of the electrostatic arcs in several regions of the phase diagram. The roots of 
the Bcthe ansatz equations for L = 64 and M = 16 are also shown. For the four smaller figures, the same 
horizontal scale is used which is shown explicitly for case (d). Cases (a) and (b) have the same vertical 
scale, as do cases (c) and (d). We note that while the roots lie close to the electrostatic arc in cases (a), 
(c) and (d), the agreement is not as good in case (b) which lies near the Moore- Read line. 



Type 


Reality 


Topology 


Parameters 


Meaning 




real 


open 


(o,e A ) eft 


Fermi pairs 




complex 


open 


e±i5 


BCS pairs 




complex 


closed 


e Al a = b 


BCS pairs 


r D 


real 


open 


(a, b) £ ft 


BEC pairs 



Table 2: Classification of the arcs containing the roots of the Bethe ansatz equation ([61]) . 



Phase 


r 


Weak coupling BCS 




Moore- Read line 




Weak pairing 


r 4 ur c u r D 


Read- Green line 


r D 


Strong pairing 





Table 3: Structure of the electrostatic arcs in the different phases of the p + ip model. 



Next we consider the solution of (I6TT) in the different regions of the phase diagram. 



4.2.2 Weak coupling BCS phase 

In this region r = T A U T B , where T A = (0, e A ), e A < to and T B is a complex arc with 
end points e ± i5. In the continuum limit, the roots belonging to T A lie between the z\. 
Recalling that the charge associated to the y m is +1, while that associated to the z\ is 
— 1/2, one finds that the effect of the arc T A can be taken into account by writing fl6Tj) as 

Jo £-y Je A e ~v y Jr B y-y 

To solve this equation we introduce a function h(y) which is analytic outside Tb, ft and 
0, and has a branch cut on the arc Tb, such that 



26 



r(y) \dy\ = —(h+(y) - h_(y)) dy, y G T B 
2m 



(64) 



where h + (y) and h-(y) are the limiting values of h(y) to the right and left of Tb- The 
appropiate ansatz for h(y) is 



h(y) = R{y) 



de 



Uo 



Q 
v 



(65) 



where 



R(y) = V(y -a)( y -b) = V(y - e) 2 + P 



(66) 



such that R(e) is positive on the interval (ua,w) and negative on (0,ua)- In particular 
R(0) < 0. Using ( 16"1|) one can write the last term of ( 1631) as 



h B y - y 



dy' h(y') 
B 2m y' -y 



(67) 



where L B is a counterclockwise contour surrounding Y B - Plugging (!65l) into (167|) and 
deforming the contour Lb around Q, and oo one finds 



dy' h(y>) 



UA ^ <p(e)\R(e)\ 

£-y 



<p{e)R{e) Q R(0) 



de 

Lb 2m y' - y J e-y J WA e-y y 

The solution of the electrostatic equation flBUl) is obtained with the choices 



I deip(e) - Q. 
Jo 



<p(e) 



\R(s) 



Q 



R(0) 



P{e) 



d£ \R(s) 



-Q 



\R(0)\ 



(68) 



Observe that the last equation in (1681) coincides, up to a rescaling of the variables, with 
the chemical potential equation ( 1541) (use |-R(0)| = \fab = \fi\). In the s-wave model the 
analogue of ( 1681) provides the gap equation. However for the p + ip model the roles of 
the gap and chemical potential equations are reversed. To find the gap equation from the 
electrostatic model one has to compute the number of roots on the arc Y: 



M 



\dy\r{y) 



\dy\r(y) 



de p(e) 



dy 
2m 



Ky) 



(69) 



Plugging ( )65|) into ( )69|) and deforming the contour Lg as done above, one finds 



ep{e) 



+ S 2 



where we have also used ( 1681) . This is the gap equation (1531) . The energy of the solution 
can be computed from (1621) 
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E 



r i 



\dy\ y r(y) + / \dy\ y r(y) = 2 dee p{e) + 



dy 



2txi 



-■ y h(y). 



(70) 



As in the case of equation (168]) . we substitute fl65|) into flTDl and after the contour defor- 
mation of Lb one finds 



E 



de p(e) e ^ 1 - 



e — e 



e-e 



S 2 



(71) 



in agreement with fl55|) . This formula is similar to the one found for the s-wave model 
except for the fact that is missing a constant term —5 2 /4 corresponding to the conden- 
sation energy [231 US]. Finally, the equation for the arc T B is given by the equipotential 
equation 



Re 



while the density of roots is 



dy' h(y' 



o, y er 



B 



r(y) = -\h(y)\, yeT B . 
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4.2.3 Moore-Read line 

This case is reached when the arc T# studied in the previous subsection closes, i.e. 5 = 0, 
so that a = b = e. The function R(y) defined in eq. (166]) becomes R(y) = y — a, so it no 
longer has a branch cut. To find the closed arc Tc one introduces an analytic function 
s(y) satisfying 



r(y) \dy\ = s(y) dy, y e T c . 
The Bethe anstaz equation is similar to ( 163]) . i.e. 



£A , P( £ ) f 
de J -^- L + / de- 

e-y Je A £-y y 



P I \dy'\4^- = 0, yyeT c 

Jvc y y 



(72) 



where is the point where Tc cuts the positive real axis. The analogue of ( 167]) is 



, »(v') 



Jt c y-y Jr c y 



u 



(73) 



where Tc is the counterclockwise contour containing the complex roots. To solve (172]) we 
choose the following ansatz 



s{y) = — 

ITT 



de 



p(e 



)_ + <to 

y y . 



(74) 
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which is motivated by the two level s-wave model where the roots form a closed loop 
when the coupling constant is less than half the distance between the two levels [231 155] - 
Introducing (1741) into (1731) one can check that (1721) holds (we use that f Tc dy'/(y' — y) = in 
for y G Tq)- The total number of roots M is given by 



M 



2 f de p(e) + [ dy s(y) = 2q 
Jo Jt c 



X — 1 



1 

9 



which shows that the criterion for the Moore-Read line is reproduced. The chemical 
potential equation is equivalent to the vanishing of the potential s(y) at the point a 
where Tc cuts the real axis, i.e. 



•, I de- 
ls e - a 



s\a) 



0. 



The latter condition follows from continuity of the weak coupling case since the density 
of roots on the arc Tb also vanishes at the end points. 
The equation of the arc Tc is given by 



Im 



dy' s{y' 



0, yeY 



c 



(75) 



where yo is a point of Tc which can be taken as a. Doing the integral one finds 



Re 



/ 



dep(e) log ( - — - 
n — 2/o 



yo 



o. 



(76) 



In ID this equation becomes, in normalized variables (i.e. y — > uy) 



Re 



VV + 1 



while in 2D it is 



- Vvolog 



+ log 



yo 



:i-Vg^ 

9 yo. 



Re 



-y log ( 1 J + y log ( 1 - — 

y) \ yo 



log 



y-i 
yo - 1 



:i-i)io g ^ 

9 yo. 



o. 



These eqs. have two real solutions, corresponding to the intersections points of T c 
with the real axis: y = a/u and y = Ea/w. The total energy can be computed in a way 
similar to eq. ([70]) and the result is 



E = 2 / dEE P (e) + I dyy s(y) = 0, 



which agrees with the mean-field result (1461) . This result can also be obtained as a limit 
of eq. (TrTj) by setting 5 = and e = a < 0. 
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4.2.4 Weak pairing phase 

This is the most complicated phase concerning the structure of the arc T which consists 
now of three pieces, T = Ta UTc UT^, where Ta = (0, ea), I\d = (a, b) and Tc is a closed 
arc that touches the points Eb = b and Ea on the real line. Combining the results of the 
two previous cases we introduce a function h(y) associated to the open arc 



Hy) = R(y) 



£-y y 



and a function s(y) associated to the closed arc Tc, i.e. 



s(y) = -R(y) 

17T 



o £ -y y 



where 



(77) 



(71 



R{y) = ^{y - o){y -b), a<b<0 
which is always positive for s G Q. The electrostatic equation (16T|) reads, 



dy> h{y') 



dE- 

o e - y Je A £-y y h D ^y' ' -y 

and it is solved by the ansatze (177] ITS]) provided 



Jl c y - \ 



y eY c or Y D 



<p{E) 



p{e) 
R(e)' 



Q 



R(oy 



and 



P{e) 



ds 1 
o R(e) 



-Q 



Q = Q' 



Qo 
R(0) 



(79) 



10) 



The last equation coincides with the chemical potential equation (154"]) since \i = R(0) = 
ab > 0. On the other hand the gap eq. (!53l) is obtained from the counting of the number 
of roots M 



I I d£ p(e) 
'o 



dy s(y) + 



dy 
2ivi 



Hy) 



ep{e) 



o ' y/(e-a)(e-b) 2G 



.(81) 



The total energy is given by 



E 



dEE P (E)+ [ dyys(y)+ [ |^ y h(y) 
'o Jv c Jl d 2ttz 

V y/(e-a)(e-b) 



(82) 



which coincides with ( 155]) . Finally the equation for the closed arc Tc is the same as (1751) 
with s(y) replaced by ([75]) . The results obtained so far must reproduce those found ealier 
for the Moore-Read line. Indeed, if a — > b, one can check that the function s(y) given in 
BJ) coincides with (J71). 
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4.2.5 Read-Green line 



This case can be obtained from the previous one taking the limit £» — ?- 0. One can easily 
prove that the number of roots on Ta and Tc is zero so that all the roots lie in the open 
arc Tp. The gap and chemical potential equations and energy are given by (j8U|8U82p 
setting b = 0. The criticality of this line can be seen from the fact that the function 
R(e)/2, which is the mean-field energy of the excitations, vanishes as \/e oc \k\ when 
| k | 0. Note that this also means that the density of roots is divergent at the origin, as 
can be deduced from (I77|) . 

4.2.6 Strong pairing phase 

In this case all the roots lie on the open arc Td = (a,b) belonging to the negative real 
axis. As in the weak pairing phase we introduce the function h(y) given by (j77j) . so that 
the electrostatic equation becomes 

p(e) q f dy' h(y>) 
/ de = / — , y G T D . 

Jo e ~y y Jl d 2m y - y 

The solution of this equation is provided by 



^ £) = WY Q = W) 



and 



Jn 



R{oy 



Notice the change of sign of the equation for Q as compared to (1751) . Computing the 
number of roots on the arc gives the gap equation: 



1 - r de <m 



2G J ' ^(E-a)(E-b) 
The total energy is given by 

dy f, / e- (a + 6)/2 



E= [ ^yh(y)= f dep{e)e(l 
Ji n 2m J \ 



y/(e -a)(e- b) 



4.3 Topological properties of the ground state 
4.3.1 Dressing and duality 

Above we established, through use of the electrostatic analogy, that the continuum limit 
of the Bethe ansatz equations for the ground state is in agreement with the mean-field 
solution. However with reference to Fig. |6l we see a noticeable discrepancy between the 
theoretical arc and the position of the roots near the Moore- Read line. Here we re-examine 



31 



the Bethe ansatz equations in the case of finite-sized systems to clarify some subtle issues 
about the nature of the ground state. Using an algebraic approach we will introduce the 
concept of dressing of states, which gives a different perspective on the duality that was 
seen in the mean-field analysis. 

For a generic splitting of the set of roots Y into nonintersecting sets Y' and Z such 
that Y = Y' U Z, we can express the Bethe ansatz equations (1231) as 



G L - L + 2M -l + Yl 



k =i y - 



E^Vm V— «■ 2|/ m . 



G- 1 - L + 2M - 1 + 



fc=l Vm 

^Vin \ *■ 2y^ 



— = — + 2^ — - — ' G z - ( 84 ) 

Setting \Y'\ = M' and \Z\ = P we take the sum over elements in Z in (1841) to give 



P (G -i_ i+2M -i) + ££^ 

2|/m V ^ 2|/^ 



2y 



Suppose that at some limiting value of G we have 

y m ^ forally m GF', 
y m = for all y m G Z. 

Taking note that M' + P = M equation (I85I) informs us that 

P = G- 1 - L + 2M => P = L - 2M' - G~ x , (86) 

while from (183]) we obtain 

CT 1 - L + 2M - 1 + £ — ^ = Y, ^^ + 2P, faey' 
^ G -i_ L + 2M'-l + ^^ 1 = £ _J^ ; y m eF'. (87) 

Equation (1H7I) is simply the set of Bethe ansatz equations for a system of M' Cooper pairs. 
The above calculations suggest that given such a solution set Y', we can augment it with 
P additional roots which all have zero value to obtain the solution set Y, provided that 
P is given by ( I86p . It is important to check that this solution set is indeed valid for there 
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are known examples where the roots of the Bethe ansatz equations must be distinct, most 
noticibly the repulsive Bose gas [36] . 

A fairly straightforward calculation, using proof by induction, leads to the following 
commutator identity for if in the form (1221) 



[if, & 



where 



PQ t C p ~ 1 (0)(GL - 1 - GP — 2GN) 



J2 Zkb k- 



k=l 



If we partition the Hilbert space of states into subspaces defined by the eigenvalues M' 
of the Cooper pair number operator N, there exists a subspace on which the commutator 
flHHj) vanishes whenever P = L — 2M' — G^ 1 which is precisely equation fl86l) . That is to 
say given any eigenstate \4>(Y')) with M' Cooper pairs, we can construct a new eigenstate 
\<p(Y)) = C p (0) \<p{Y)) of M = M' + P Cooper pairs and the same energy as \(f>(Y')) 
provided ()86[) holds. Whenever this is the case we say that the state \(p(Y)) is a dressing 
of the state \(j)(Y')). 

To explore in more depth the consequences of this result, we first note that P must 
necessarily be a non-negative integer. In terms of the filling fraction x' = M'/L, non- 
negativity of P implies that 

2 V 9. 

so that the state \<j)(Y')) belongs to the strong pairing phase. Letting x = M/L 
(M' + P) /L denote the filling of the dressed state we find 



x + x 



1 - 



9 

which is simply the duality relation ( 1501) obtained from the mean-field theory. We can 
immediately conclude that \<j)(Y)) belongs to the weak pairing phase. However one signif- 
icant point about this derivation of the duality is that it only applies in particular cases 
when the coupling g is rational. This aspect has some ramifications to which we will 
return later. 

For the moment, we will consider a simple example of dressing. Starting with the 
sector of the Hilbert space where the number of Cooper pairs is zero there is only one 
state, that being the vacuum |0). We can dress this state with M additional Cooper pairs 
where, by equation (1861) . we must choose 



M = P = L — G~ 



39) 



There is no reason yet to expect that the dressed state will be the ground state in the 
sector with M Cooper pairs. However we do recognise that (1891) is precisely the equation 
of the MR line f !45p found from mean-field theory, for which the mean-field ground-state 
energy is zero. The dressed state will have the same energy as the vacuum, which is also 
zero. Performing the dressing gives us 



C M (0)|0> 



Eh 

,k=\ K 



n M 



keK+ 



°k L -k 



-I M 



10) 
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which is precisely the MR state (l4"3"j) . Recall that the projected mean-field wavefunction 
gives the MR state as the ground state on the MR line with zero energy. We can show 
that this is indeed the ground-state. In fact, in all instances dressing a ground state in 
the strong pairing phase leads to a ground state in the weak pairing phase, which we now 
show follows from the Perron- Frobenius theorem. 

By adding an appropriate constant term to (1221) we obtain an operator with strictly 
negative entries for each sector of M' pairs. The Perron-Frobenius theorem tells us that 
the ground-state vector of this operator can be normalised such that all components of 
the vector are strictly positive. Applying the dressing operator C p (0) with P given by 
( 1861) to the ground-state of M' pairs leads to an eigenstate of ( 1221) in the weak pairing 
phase, which also has strictly positive components. By the Perron-Frobenius theorem, 
this must be the ground state in the weak pairing phase for M' + P pairs. 

It is of importance to compare this result with the previous discussions of the solution 
in the thermodynamic limit described in Subsection 14.21 There, it was concluded that 
on the MR line the ground-state roots formed a closed arc in the complex plane. In the 
finite system, we have now found that all ground-state roots are zero on the MR line. In 
both instances the ground-state energy was found to be zero. One way to resolve this 
apparent contradiction would be to suggest that the assumption used in Subsection I4.2[ 
viz. that in the thermodynamic limit the ground-state roots form a dense arc in complex 
plane, was not valid. To investigate this matter further we next turn our attentions to 
the numerical solution of the Bethe anastz equations. 

4.3.2 Numerical solution of the Bethe ansatz equations for finite systems 

First we will discuss the behaviour of the roots in the three different regions. Later 
we will discuss properties across all regions. The numerical solutions of the Bethe ansatz 
equations are obtained starting from the initial condition y m — > (1 + G)k 2 (m = 1, . . . , M) 
as G — > 0, with the k chosen to fill the Fermi sea. As g increases, the roots y m , closest to 
the Fermi level become complex pairs. When g approaches the MR-line the roots bend 
towards the origin, and at the value x = xmr all the roots collapse onto the origin. At 
larger values of g one enters the weak pairing phase where all the roots are non-zero, 
except at some rational values of g where a fraction of the roots collapse again. Finally, 
in the strong pairing regime all the roots become real and they belong to the interval 
(a, b) on the negative real axis. The RG-line is obtained when b = 0, in which case the 
quasiparticle energy f l3"§j) becomes gapless. 

In the weak coupling regime the parameters a, b appearing in (|47|) give the end points 
of the complex arcs shown in Figs. UJa) -[7(d) for the ID and 2D models (recall Table 1). 
On the MR line a and b become real and equal, giving E = from (1481) . and the complex 
arc closes satisfying eq. (ITS]) . 

Fig. [TJe) shows how the roots approach this closed arc for increasing values of M at 
^ = 1.99 an d x = 1/2. For g = 2 and x = M/L = 1/2 (MR line) all roots collapse to the 
origin in any finite system. The continuous arc is the solution of eq. f J76|) corresponding 
to g = 2,x = 1/2, where L — > oo is taken before g — > 2. However the exact solution 
of the Bethe ansatz equations f[2"§|) is y m = 0, Vm, exactly giving the Moore-Read state 
[4"5] C(0) M |0) which also has total energy E — 0. In Fig. [7(f) we plot the real part 
of the Bethe ansatz equations roots for a 2D momentum distribution in a system with 
M = 8, L = 32. The Moore-Read points, in which all the roots collapse to zero, are 
clearly noticeable. This fact presents a paradoxical situation where the thermodynamic 
limit does not commute with the limit where the coupling g approaches the MR line from 
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the weak coupling BCS phase. 

In the strong pairing phase we expect from the calculations in the continuum limit 
that the ground-state roots should be real and negative. Moreover, the Perron-Frobenius 
theorem indicates that there can only be one such set of roots with this property. This 
fact makes the numerical analysis of the ground-state roots fairly efficient in the strong 
pairing phase, since the solution converges quickly to the unique real and negative solution 
set. In all instances we considered we always found a numerical solution set of real and 
negative roots in this phase. 

To illustrate the results for this phase we define the discrete root density to be 

(M — — Uj) 

which we have plotted in Fig. |8]for a system with L = 450 and M = 150. The results 
compare the cases g = 3, the Read-Green line, and g = 20. As discussed in Subsection 
14.2.61 the continuum limit of the Bethe ansatz equations predicts a divergence in the root 
density at the origin when the coupling is on the Read-Green line. The numerical results 
support this picture with the value of the root density at the origin for g = 3 several times 
larger than for g = 20. 

In Fig. Owe plot the exact and mean-field energy of the ground state of a 2D system, 
showing results for the two filling fractions x = 1/2, 1/4. Observe that at the MR coupling 
the exact and mean-field energies coincide, i.e. E = 0. For x — 1/4 there is a departure 
between the exact and mean-field results in the strong pairing phase, which disappears 
for large systems. 

In Fig. ITUT a) we plot the expectation values of the Cooper pair density (N n ) (using 
11171 in Appendix I A. 2 1) . In Fig. ITOT b) we plot the exact and mean-field values of the the 
fluctuation parameter C defined as 



In both of these instances there is reasonably good agreement between the mean-field and 
the exact results. 
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Figure 7: Real and imaginary parts of the roots of the Bethe ansatz equations in several regimes for 
several values of the coupling constant g. In the cases (c,d,e) the corresponding electrostatic arcs are also 
shown: (a) ID model, M = 2,1,1 = 62, (b) ID model, M = 16, £ = 64, (c) 2D model, M = ?A.L = 62, 
(d) 2D model, M = 16, L = 64 (e) 2D model, g = 1.99, M = 8, 14, 31, L = 2M. Notice that the MR point 
is g = 2, x = 1/2. (f) Real part of the Bethe ansatz equation roots for 2D momentum distribution with 
AI = 8, L = 32. The range of g covers both the weak pairing and strong pairing phases. The dressing 
points, in which some of the roots collapse to zero, can be seen to occur at regular intervals between 
9mr =4/3 and g RG = 2. 



36 



rfy) 

0.000035 ■ 
0.00003 ■ 

0.000025 ■ 
0.00002 ■ 

0.000015 1 



• 


g=10.0 




g=15.0 




g=20.0 



0.00001 ' 

-14000-12000-10000 -8000 -6000 -4000 -2000 



r(y) 

0.012 r 



0.010 
0.008 
0.006 
0.004 
0.002 
0.000 



-600 -500 -400 -300 



-200 



-100 



Figure 8: Discrete root densities as given by ([90]) for L = 450 and M = 150. The momenta have been 
chosen as |k| = m, m = 1, ...,450. The left graph shows results for the coupling choices g = 10, 15,20, 
which correspond to the strong pairing phase. The support of the discrete root density lies on the negative 
real axis, and in each case the maximum value does not occur at the endpoint. For comparison, the right 
graph shows analogous results for g = 3 corresponding to the Read-Green line. Here the maximum value 
of the discrete root density occurs at the endpoint of the support at y = 0. Note also the different scale 
for the vertical axis. 




Figure 9: Comparison between the exact and mean- field ground state energy in 2D for x = 1/2, M = 
14, L = 28 (left) and x = 1/4, M = 16, L = 64 (right). The arrows indicate, from left to right, the values 
of the coupling g corresponding to the Volovik, Moore- Read and Read-Green lines. For x = 1/2 the latter 
point is absent. 



4.3.3 A zeroth-order quantum phase transition 

We have already observed that the MR line has the peculiar property that the behaviour 
of the roots of the Bethe ansatz equations in the limits g — > Qmr-, L — > oo depends 
on the order in which the limits are taken. Another curious property of the MR line is 
the discontinuity of the ground-state energy E(g,x) in the thermodynamic limit as the 
filling fraction x approaches the value xmr from the weak pairing region. To derive this 
result, for finite L we take the one-pair state and dress it to give the dual GS in the weak 
pairing region. The filling xj of the dressed state is given by (|50|) . setting xjj = 1/L, 
i.e. xj = xmr — which approaches xmr = 1 — l/g as L — )• oo. Since the MR pairs 
carry no energy, the ground-state energy of the dressed state coincides with the one-pair 
energy. To compute this energy we consider the Bethe ansatz equation for one Cooper 
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Figure 10: a) Cooper pair expectation values of the n-th level, (N n ), for a 2D momentum distribution 
system with M = 24, L = 48, for several values of g. The values of g corresponding to the weak coupling 
BCS phase. The mean-field counterpart of this magnitude is also shown for each value of g and x=l/2 
(continuous lines), (b) The parameter C as given by (p?Tj) . measuring fluctuations of the expected value 
of the ocupation number is plotted as a function of g. Results are shown using the roots of the Bcthc 
ansatz equations for a M=8, L—32 system, and for mean-field case in the continuum limit with x = 1/4. 
The arrows indicate the Volovik, Moore-Read and Read-Green points. 



pair and take the continuum limit (i.e. eq. ( 121?]) with M — 1). For simplicity we take the 
momentum distribution to be that for free particles in 2D, i.e. p(e) = L/oj. This leads to 

T i , y 1 f u de 1 

L- — -l= > — =>l-- = y / . 

G y - k 2 g J u y-e 

This equation has a unique negative energy solution y < satisfying 

log 



g u \y-uj 

which we denote as y = S(g). From here one derives the aforementioned discontinuity on 
the MR line xmr = 1 — <7 _1 , 

lim E{g,xj) = £{g) ^ E(g,x MR ) = 

L— >oo 

which may reasonably be termed a zeroth-order quantum phase transition. We note 
that in thermodynamic systems analogous zeroth-order transitions have been discussed 

in pzgCESlHE]. 

All of the above considerations point towards the MR line signifying a phase boundary 
which has significantly different properties to the RG line. In the latter case, both the 
excitation spectrum is gapless and the mean- field fidelity susceptibility is divergent. These 
are not properties associated with the MR line. 

4.3.4 Winding number of wavefunctions 

Here we formulate a description of the three phases in terms of the ground-state wave- 
function topology. To this end we will use a winding number approach. 
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Let us consider a continuous map S 2 — > S 2 , given by a function g : (X, Y) — > (g x , Qy) 
where (X,Y) and {qx,Qy) are the stereographic coordinates of the spheres. The winding 
number is defined as 



w 



dQx dg Y 



The function 



(1 

dJ($x,QY 
d(X, Y) 



Qx 



Qy) 2 



dX dY 



9J(qx,Qy] 



d(X,Y) + + 



(92) 



is the Jacobian of g. An alternative way to write (1921) is 



w 



-I dX dY 



dgdg _ dgdg 

dz dz dz dz 



(93) 



where g = gx + i&Y and z = X + iY . 

To start, we consider the exact eigenstate of ( JT9l) for one Cooper pair, as given by ( |28l) . 
In the limit L — > oo with one Cooper pair there are only two ground-state phases, the 
weak coupling BCS phase and the strong pairing phase, with the transition point given 
by g — 1 at which the ground-state energy is zero. In first quantized formulation we have 
that the momentum-space wavefunction is 



0(k) 



k 2 -E 

where E denotes the energy. In terms of the complex variable z = k x + ik y , the Jacobian 
is given by 

dg dg dg dg (zz) 2 — EE 

dz dz dz dz 
The winding number is then obtained as 



J 



w 



du- 



[zz - E) 2 (zz - E) 2 

(zz) 2 - EE 
[zz + (zz — E)(zz — E)} 2 
u 2 -EE 



which yields 



u 



w 



[u + (u - E)(u - E)} 2 
u 

(E + E)u + u 2 + EE 



E = 0, 
E^O. 



(94) 



In this manner we find a topological change in the wavefunction when the energy is zero, 
or equivalently g — 1, consistent with the phase diagram Fig. [2j However in both the weak 
coupling BCS phase and the strong pairing phase the wavefunction topology is trivial. 

It is useful to contrast the above with the analogous result for the mean-field wave- 
function. Up to a factor of A* which can be omitted, we have from (I41|47l) 



g(z,z) 



e(z, z) — zz + n 



i(z, z) = \J (zz — a)(zz — b) 



where fi 2 = ab. The Jacobian is given by 
1 



J 



(2zz e(z,z)f 



[((a + b)zz - 2ab - 2/2 e(z, z)) - (zz) 2 (2zz -a-b- 2e(z, z)) ] 
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We have found numerically that 



1, /i>0, 

w = <( a /(a-l), /i = 0, (95) 
0, fi < 

where = corresponds to a choice a < 0, b = 0. In this case the winding number is not 
an integer at ji = . The reason for this is that setting ji = gives 



g(z, z) = (Vii - a - \J~zz) 



which is not continuous at z = 0. The above shows very different topological predictions 
between the exact wavefunction and the mean-field wavefunction, specifically with the 
mean-field wavefunction having non-trivial topology in the weak coupling BCS phase. 

To gain a further understanding of the many-body system, next we consider the exact 
wavefunction for two Cooper pairs: 

,/,/, = „ = \ _/ Q(ZuZ 1 ,E 1 )q(z2,Z2 j E 2 ) +B(Zi,Zi,E 2 )g(Z2,Z 2 ,Ei), Z!^Z 2 , , QR x 
V\Z\, Zi, z 2 , z 2 ) — < q Zi = Z 2 

where g(z,z,E) is the one-pair wavefunction and Ei, E 2 are roots of the Bethe ansatz 
equations ( 15U|) . which can either be real or a complex conjugate pair. In order to apply the 
winding number approach described above we reduce f )96|) to a complex-valued function 
by setting 

q(z, E u E 2 , c) = g(z, E 2 , E x , c) = ip(z, z + c, E u E 2 ) 
for some 0. Numerically we find 



2, E x = 0, E 2 = 0, 
w={ 1, E x = 0, E 2 ^ 0, 
0, Ex^O, E 2 ^0. 

This result shows that the winding number coincides with the number of MR pairs. 
Equation f[9Uj) is given an interesting interpretation if we parameterize g(z, z) as in the 
BCS model 



g(z, z) = V ^, Z,Z }. , u(z, z)u(z, z) + v(z, z)v(z, z) = 1 
where u and v are well defined functions in the complex plane. It can be shown that 

w = — ^— / dX dY [d x {u d Y u + v d Y v) — dy{u d x u + v d x v)] . (97) 
and from the Gauss theorem 

/■2tt 

w = lim / - — : (u dgu + v dgv) (98) 



r— >oo 



2tu 



where X + = re* 9 . This expression indicates that it) measures the angular dependence 
of the wave funcion g(z,z). The result ( |95l) for the BCS wavefunction can be derived 
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analytically from (I98p . by noticing that the proper definition of u and v depend on the 
sign of the chemical potential, i.e. 




zA 

u = _ _ _ , v =\l — — 3; j /x > 0, 

\/2e(z, z)(e(z, z) - zz + /x) 

le(z,z) + zz- fi 

U = \ ; — , V = , /X < 

V 2e{z,z) y/2e(z,z){e{z, z) + zz - /x) 

which leads to 

ty = lim \u(r, 0)| 2 = 1, /x > 0, w = — lim \v(r, 0)| 2 = 0, /x < 0. 

j — >oo r—>oo 

Similarly, equation (194)) can be derived from ( 198]) using the following expression of the u 
and ?; functions corresponding to the one pair wavefunction 

z 1 
u = v = E = 0, 

y/l + ZZ \/l + ZZ 

zz — E z . 

u = _ _ , v = _ _ , E ^0. 

\J\zz — E\ 2 + zz ^J\zz — E\ 2 + zz 

In the general case of an exact wavefunction with M pairs, having P vanishing rapidities, 
the asymptotic expression of u and v immediately yields the expected value of w, i.e. 

u ~ z M z M ~ p , v~z M ~ p =^ W = P. 

This result shows that non-trivial topology of the wavefunction can only be found in the 
weak pairing phase. 

As a final comment we note that the winding number w appears in a completely 
different context, namely the 0(3) non-linear sigma model in 2D. The order parameter of 
the latter model is a three component unit vector n, which can be constructed out of the 
two component spinor ip = (u,v), as n = ip^aip. Using this parametrization the winding 
number w can be written as 

w — — / <i 2 x n ■ {dJi x d u n) e^ v 
8tt J R 2 

where e^ v is the 2D Levi-Civita tensor. This term appears in the action of the 0(3) non- 
linear sigma model multiplied by a parameter 9, which is defined modulo 2n. In a famous 
work, Haldane mapped the ID antiferromagnetic Heisenberg model of spin S into the 0(3) 
non-linear sigma model with 9 = 2nS [30] ■ The gapped nature of the integer spin chains 
then follows from the gapped nature of the non-linear sigma model with 9 = (mod 2ir). 



4.4 Vortex wavef unctions 

So far we have discussed the p + ip model for a constant value of the order parameter 
A, which describes the system in the bulk. It is however of interest to study the effect 
of vortices. As shown by Read and Green [53] there are vortices with zero energy when 
the chemical potential /x is positive. We shall show below that the structure of the 
vortex wavefunction for /x > depends on the region of the phase diagram, with different 
behaviours in the weak coupling BCS and weak pairing regimes. However the transition 
between these two regimes is smooth in the mean-field analysis. 
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We shall consider the Bogolioubov-de Gennes (BdG) equations in the simplest situa- 
tion of a single vortex with rotational invariance characterized by an order parameter 

A(r) = ie** |A(r)| 

where I is the vorticity and r, tp are the polar coordinates centered on the vortex. (We 
shall closely follow the results of references [291 EH E9] with the appropiate notational 
changes: the mass is equal to 1, the gap A and the chemical potential in this paper are 
twice those considered in the cited references.) The BdG equations are a set of coupled 
equations for the functions u n (r), v n (r), corresponding to energies E n , which determine 
the Bogolioubov quasiparticle operators j n by 



In 



where a(r),a'(r) are annihilation and creation operators of polarised electrons. The zero 
energy solutions, i.e. E n = 0, can be choosen to satisfy the reality condition u n (r) = v n (r), 
in which case the quasiparticle operator 7„ corresponds to a Majorana fermion. For odd 
vorticity, I = 2n — 1, the BdG equation for u(r) can be most conveniently written in terms 
of a function x( r ) defined through 



u(r) = x(r) exp (~ J dr' |A(r')M 



and it reads 

„ x' /|A 2 (r)| n 2 



Let us suppose that the order parameter |A(r)| vanishes inside a core of radius a and that 
it reaches its bulk value, A , outside i.e. 



|A(r) 



0, < r < a, 
A , a < r < oo. 



The previous equation becomes 



r x" + r x' + 



2 X" + rx + [^r 2 -n 2 ]x =0, < r < a, (99) 
AST 



»-^)r 2 -n 2 



X =0, a < r < oo. 



Notice that the coefficient multiplying the r 2 term in the second equation changes sign 
when crossing the MR line fi = Aq/4. This fact leads to two different type of equations 
corresponding to the weak coupling and weak pairing regimes. (We do not consider the 
strong coupling regime where the vortices have a trivial topological structure [2H1 1531 159"].) 

4.4.1 Weak coupling BCS phase 

For the weak coupling BCS phase fi > Aq/4 we make the change of variables 

A 2 

,, — if 2 a - — ° - k 2 
r — A o> r 4 — ' 

Eq. (1991) becomes 

x 2 Xxx + xxx + [x 2 -n 2 ]x =0, x = k r, < r < a 
y 2 X yy + VXy + [y 2 - n 2 } X =0, y = kr, a < r < oo 
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whose solutions are Bessel functions. Let us denote by I and II the regions < r < a 
and a < r < oo respectively. Imposing that x is regular at r = and infinity one has 

Xi = AJ n (k r), 
Xii = BJ n {kr) + CY n {kr). 

The constants A, B, C are found by matching the function u(r) and its derivative u'(r) 
at the boundary r = a. The relation between u and x m t ne regions I and II is 

ui(r) = Xi(r), 

u n {r) = X ii(r)e- Aor/2 . (100) 

Hence 

u I (a) = u H (a) xi{ a ) = X//(«)e~ aAo/2 , 
<(a) = < x (a) =► X '/(«) = (XW(«) " ^X//(a))e- aAo/2 - 
Combining these two equations one can write 

XM) = xiM)e- a ^l\ 
X» + ^Xi(a) = x'ii(a)e- aA ^ 2 

with the explicit expressions 



AJ n (k Q a) 
k J' n (k a) + ^-J n (k a] 



(BJ n (ka) + CY n {ka))e- aAo/2 , 
k{BJ' n {ka) + CY;{ka))e- aAo/2 . 



This then leads to 
A 2 



-a7re aA °/ 2 



kY^(ka)J n (k a) - Y n (ka) ( k J' n (k a) + ^ J„(fc a] 

k J' n (k a) + ^J„(/c a) 



For an infinitesimal core, a — > 0, the zero mode is localized on an isolated vortex with 
wavef unction 

u (r) = J n (kr) e" A ° r/2 , r > (// > A 2 /4). 

This result reproduces the one found by Gurarie and Radzihovsky [22]. The funcion u(r) 
has an oscillatory and decaying asymptotic behaviour dominated by A , i.e. 



u\r) 



cosjkr + O ) _ rAo/2 



r — > oo. 



4.4.2 Weak pairing phase 

To study the weak pairing case < /i < Aq/4 we make the change of variables 



A 2 
4 
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which corresponds to the replacement k — y ik in the previous case. The equation for x 
becomes 



x 2 Xxx + xxx + [x 2 - n 2 ] x =0, x = k r, < r < a, 
V 2 Xyy + VXy ~ [y 2 + n 2 ] x =0, y = Ax, a<r < oo. 

whose solutions are Bessel functions in region I and modified Bessel functions in region 
II: 

Xi = AJ n (k r), 
Xii = BI n (kr) + CK n (kr). 

The computation of A, B, C is similar to the previous one. The final solution is 



B 
A 
C 
A 



aA /2 



-ae 



_ ae aAo/2 



kK' n (ka)J n (k a) - K n (ka) ( k J' n (k a) + — J n (k a) 
kl' n (ka)J n (k a) + I n (ka) ( k J' n (k a) + ^ J n (k a) 



(101) 



In the limit a — > one obtains 

u{r) = I n {kr) e - Aor/2 , r > (fx < A 2 /A) 

in agreement with [29]. It is important to notice that the exponential factor e -A ° r / 2 in 
equation fllOOp is essential to get a convergent asymptotic behaviour of u(r). Assuming 
that the coefficient B, given in eq. fllOip . is non zero one finds for B ^ 



u(r) 



1 



~ — =exp 

r 



A 2 



r — y oo, 



while if B = 0, the asymptotic behaviour of u(r) changes to 

/ A 



u\r) ~ 



exp 



A 2 
_£ 

4 



r —y oo. 



which is a more localized vortex. The condition for this to happen is 

B = 



K n (ka) 



J n (k a) 



Ao 
2 ' 



(102) 



This equation has an infinite number of solutions for a as a function of \i and Ao- In the 
case where ka, k$a » 1, the solutions of eq. (11021) are 



i mr fk + A /2 

kna ~ — I h arctan 

4 2 \ k 



(mod 7r). 



Hence there is a minimal core size in order to have these more localized vortices. In Fig. 
[TTJwe plot the function u(r) for several values of the coupling constant g. In the weak 
coupling regime the function u(r) has an oscillatory behaviour, while in the weak pairing 
regime the oscillations are absent. We observe that the different qualitative behaviour 
depends on the region of the phase diagram. 
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Figure 11: The wavefunction u(r) for the vortex solution with I = — 1 for £ = 1/4 and several values 
of the coupling constant g. The size of the core is a = 6 in units of ^fuj. Left: weak coupling BCS region 
for g = 0.3 — 0.13 in steps of 0.2. Right: weak pairing region for g = 1.4 — 1.9 in steps of 0.1. 

5 Conclusions 

In summary, we have constructed an exactly solvable pairing model which contains in 
various limits some of the most studied pairing models of current studies such as: the 
Richardson model with s-wave symmetry, the p + ip BCS model, the Russian doll BCS 
model and the Gaudin models with s-wave and p-wave symmetries. We have computed 
the wavefunctions scalar products and one-point and two-point correlators. Futhermore, 
we have carried out an in depth study of the p + ip model combining the exact solution 
with mean-field methods. We have found the phase diagram of this model characteriz- 
ing the various phases in terms of the fidelity susceptibility, topological invariants, and 
other quantities such as the ground state energy, gap, chemical potential and occupation 
numbers. We have also studied the structure of the vortex wave function which is shown 
to depend on the region of the phase diagram. The richness of the p + ip model calls 
for further studies concerning its behaviour. This task can be greaty facilitated by the 
existence of an exact solution for studies such as entanglement [17], correlation functions 
[31 HE] and quenching (BJ F20]. 
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A Calculations for the anyonic pairing model 

A.l Cyclic renormalisation group 

In this Appendix we discuss how the Hamiltonian (jSJ) is related to a known Hamiltonian, 
viz. the Galzek- Wilson model [2U [25]. This model was introduced as a simple example 
of a quantum mechanical system which admits a cyclic renormalisation group (RG) map. 
As will be shown below, the Glazek- Wilson is simply the one-body version of the anyonic 
pairing Hamiltonian flSJ). Cyclic behaviour of the RG also occurs in the many-body system, 
and we will demonstrate how the RG equation is easily obtained from the Bethe ansatz 
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equations. 

The model introduced in [!Ml [25] . which we will denote as Hqw-, is defined by the 
matrix 

H GW (n, m) = \j£ n £ m {6 n)m - g - ih sign(n - m)) (103) 

where M. < n, m < Af. In [2U [25] the £ n are assumed to scale exponentially £ n = /i n , \i > 
1 under which there is a discrete scale invariance n — > n + 1 such that 

H{n + 1, m + 1) = /i m). 

However this scale invariance is broken by the infra-red and ultra-violet cutoff's A4 and 
Af. 

In [2U [2S] the RG is performed by a Gaussian elimination procedure for the highest 
component of the wavefunction in the Schrodinger equation, reducing Af to Af — 1. Under 
this process the result of the RG is that the coupling h remains invariant while the coupling 
g renormalises as 

9N-i = -: • 104 

1 QN 

Reparameterising the couplings in terms of angles 0^, x through 

tan(#/y) = 

h 

tan(x) = h 

the RG equation becomes 

, , , s tan(0 Ar ) + tan(x) , , , 

tan(0^-_i) = — 7—^1 — 7-r = tan(0^ + xj 

1 — tan(0_^) tan(x) 

or equivalently, iterating the RG equation p times, 

<t>M-p = <\>N + PX- 

A cycle occurs in the RG after p steps if \ = n/p since 

9Af- P = /itan(0 A r_ p ) = h tan(0 A r + px) = hta,n((j) A f + 7r) = g^ 
with the RG period given by 

7T 

P 



arctan(/i) 

Next we show that (11031) is equivalent to the one-pair version of (JS]). The first thing 
to observe is that for the one-pair version of (jSj) the string operators Q appearing in the 
definition of the anyonic creation and annihilation operators can be removed by a unitary 
transformation. The Hamiltonian is then given in terms of the hard-core boson operators 

as 

^one-pair = ^ ZjNj ~ • 7~^m X] ^ r ( eX P(-^ a )^ & fc + h.C.) 
0=1 v y k>r 
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with the Hilbert space spanned by the one-body states 

\j) = b]\0), j = !,..,£. 



Here the energy is 



Eone— pair ri \ 2/ 



sm(a) 2 
sin(a — 2/3) 



such that 



exp(-2za) J] / "YlK = 1- 
J -- L 1 - o l zi 



Setting 



yields 



£ k = z 2 k , g = tan(2/3) cot (a), /i = tan(2/3) (105) 



= (1 - #)#one-pair, N - M = L. 

To derive the RG equations from the Bethe ansatz, we only assume that the variables 
Zf. are distinct and ordered < %\ < < ... < Zl- We may then write 



exp(-2m) 1 ~ q2y TT l ~ < i 2 y 2z l = l 
' 1 - q~ 2 y~ 2 z 2 L " Al 1 - q- 2 y~ 2 zl 



Assuming >> which is a reasonable assumption for low energy states in a weakly 
coupled system, we can eliminate the highest energy degree of freedom of the Hamiltonian 
by making the approximation 



1 - q'y 2 z 2 L ^ 4 



1-q 2 y 2 z\ 

exp(-2za + 4z/3) TT - - 9 | w 1. 

AA 1 — a~ 2 v~ 2 zr 



q 'y >z% 

Thus for one RG step we have the renormalised coupling 

a/V-i = QtjV — 2/3, (106) 

while /3 remains invariant. Substituting back into fl 1 5 [) shows that, in terms of the Hqw 
coupling parameters, h is invariant while 

9M-1 = tan(2/3) cot(a^_i) 

= tan(2/3)cot(a^ - 2/3) 

cot(ajv') + tan(2/3) 



tan(2/3)- i _ tan ^ cot ^ a ^ 

9N + h 2 
1 - 9M 



in agreement with ( 1104ft . 
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For the many-body Hamiltonian (jSJ), a completely analogous approach can be taken 
to establish the existence of the cyclic RG equation. Starting with the Bethe ansatz 
equations (llip we can make the same approximation as in the one-pair case by eliminating 
the highest energy level, i.e. if \zl\ » \y m \ for all m = 1, ...,M then 



1 - g l y n }z\ 



Vm 



M 



and ffTTj) gives 

exp(-2za + 4i(3) ] 



L-l 



k=l 



i - q 2 yJ4 
q~ 2 ym z l 



M 



1 - 9 4 2/ m 2 2/f 



Q~ 4 ym 2 y] 



m 



1,...,M. 



This leads to exactly the same RG equation (11061) for a, with /3 again invariant. 

A. 2 Wavefunction scalar product and correlation functions 

One of the byproducts of formulating the model via the QISM is that it readily provides 
access to the evaluation of correlation functions. Again, we will only outline these calcu- 
lations. We will follow the general procedures of [ID] , where explicit details of derivations 
can be found. 

It is convenient to first rescale the L-operator as 



/ 1 



L(x) 



such that it has the properties 











q q ^x-qx 

qx—q~ 1 x~ 1 



q x—qx 



\ 



\ 




1 / 



(107) 



L(x] 
L(x] 



9=1 



x=q 



I® I 
P 



with P the permutation operator. We also rescale the matrix U to be 

~ _ / exp(— 2ia') 
\ 1 

These rescalings have the net effect of working with rescaled elements of the monodromy 
matrix 



T](x) 



exp(— ia 



n 



^ qxz k 1 - q 1 x 1 z k 



The transfer matrix eigenvalues of t(x) = T\ (x) + T 2 2 (x) are now 
A(x) 



M 2 -1 -2 -1 

~ / \ 1 r q xy j - q yjx 
a(x)[[ 



xy 3 



i 



y s x- 



-i 



A/ -2 -1 2 -1 

d{x))Yl 



\ %yj 1 - Uj-x 1 
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with 

L -1 -i _ -l 
~f \ i o- i\ TT " XZ k " x Z k 
a(x) = exp(-2ta ) I I ZI — — 

1 xz k ~ % x z k 

d(x) = 1 

One reason for working with the the rescaled monodromy matrix is that it provides a 
compact solution of the "inverse problem" of expressing the local operators bj, Q, Nj in 
terms of the elements of the monodromy matrix. The result, which is obtained by a minor 
generalisation of the methods used in [HJ EOJ H7] , is 

h = ( n f (9*o ) ^fe) (n f ~v* 

h ) = { Il^Zk) J Tl{q Zj ) (f[r\qz k 

\k=l J \k=l 

i-nj= (n%**) ) 2?(9%) ( n* _i (?2* 
\fc=i / \fe=i 

Note that the operators {t(qzj) : j = 1, L} are conserved. Explicitly these are 

t(qz k ) = L k{ k-i)(qz k zj:\)...L kl (qz k Zi 1 )U k L kL (qz k z^ (108) 

where 

LkiiqZkZf 1 ) = qZkZ^ + q 1 z k 1 z l /2N k +2N t -2 ^ g 2-2jv fc -2jv ; \ 

+ 2 g i~ g " 2 2 -i (blbi + b k b\ 
q l z k z x -q 2 z k zi V 



(e - q~ 1 )(q 2 Zkzr 1 - q^z^zi) 



The eigenvalues of t(qz k ) are 



M _i -1 -1 

A (^fc) = I rr T =T- ( 109 ) 

f^QZkVj -q~ l yjz k 
Next we consider the scalar products between the states (<&(W)|, |$(F)), where 

L 

(*(w)\ = (oin^K), (no) 

L 

my)) = U f ^\°y ( in ) 



If iy = {u>j} is a set of solutions of the Bethe ansatz equations ffTTj) . and F = {^} are 
arbitrary parameters, then the Slavnov formula for the the scalar product is (57j 



($(W)|$(F)) =C(W,Y)det(F(W,Y)) (112) 
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where elements of the matrix F(W, Y) are given by 



F(W,Y)i 



(q 2 - q 2 )d(wj) 
{VjWi 1 - y7 l Wi) 



M 



hi 



and 



C(W,Y) 



+d(yj) n^^v 



1 



Vj 1 w m q 2 ^ 



;ii3) 



Note that 



C{Y,Y) 



Uk>i(ykyi 1 - v k 1 yi)( w i w k 1 - w i V) 
i 



(114) 



n^li U^ P (y P y q 1 - y P %) 

We remark that the Slavnov formula (1112jl is a scalar product as opposed to an inner 
product. However for the anyonic pairing model we expect, due to a lack of symmetries, 
that generically there are no degeneracies in the energy spectrum for fixed particle number. 
In such an instance the eigenstates must be real, up to an overall phase factor, in which 
case the scalar product between any two eigenstates is equivalent to the inner product up 
to a phase. 

An application of the Slavnov formula (I112p is that it allows us to compute the nor- 
malised wavefunction overlap between two states with different values of the coupling 
parameter a. Specifically this permits us to compute the fidelity, as defined in [651 168] - 
which has been proposed as a technique to identify quantum phase transitions (for a 
review see [2S]). To make this point transparent, we start with the explicit form of the 
monodromy matrix elements in terms of the elements of the L-operators which is 



■nix) 



E 



9 j 



L 1 1 L ~ 1 \ 



k 1: ...,k L =l 

Since g is a diagonal matrix we can then write 



E 

fel,...,/C£_l— 1 



and each of the operators 



Lkl( XZ 2 \ 



J J 



[XZ-, 



(115) 



(116) 



is not explicitly dependent on a. Consequently the explicit a-dependence of Tj(x) is 
only in terms of an overall phase. Thus we can take the set W to be a solution of the 
Bethe ansatz equations for a given coupling, say ai, and the set Y, which in general 
may be arbitrary, to be a solution to the Bethe anstaz equations with a second choice 
of the coupling, say a 2 - This is not to say that the eigenstates of the Hamiltonian (JSJ) 
are independent of a, but rather the dependence is implicit through the solution sets W 
and Y of the Bethe ansatz equations. In this way the normalised wavefunction overlaps 
between ground states with different values of the coupling a can be computed, which is 
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in essence the computation of the fidelity^ Note however we cannot compute the overlaps 
between ground states with different values of q nor the Zk, since these parameters appear 
explicitly in the operator elements (I116p . 

Using the Slavnov formula and the solution of the inverse problem the form factors 
(matrix elements) for the local operators bj, bj can also be written, although we will not 
give details here (however see Appendix IB.4I for an example). Our main objects of interest 
for now are the one-point correlation functions, or expectation values, of the Cooper pair 
number operators Nj . We begin with the following form factor between two Bethe states 

m 

($(w)|T|(*)i$on> 

M / -1-2 -1 2\ 

= C(T^y)d(x)0(z)TT^T-^ X i V det(F(W,Y) + Q(W,y;x)), 

where 

fc = i x V k ~ x Vk 
and Q(W, Y; x) is a rank-one matrix with elements 

(q 2 - q- 2 )d{wi)d{yj) 



Q(W,Y;x) 



(xw i 1 — x 1 w i )(xw i x q 1 — x 1 w i q) 

M / ~ / s ivi 

n1 9 | 9\ / cl(x) -r-r xw, q — X Wuq 



This expression enables us to compute the expectation value 

($(Y)\N m \Q(y)) 



(N m ) 



>MY)MY)) 

(<b(Y)\f 2 \qz m )mY)) 



1 - t 1 (qz r 



(*oo\*<y)) 



det(F(Y,Y) + Q(Y,Y ]q z m )) 

det F{Y,Y) " 1 ' 

Note that for the present case because d(qz m ) = 0, d(x) = 1 we have the following 
simplification for the matrix Q(Y, Y; qz m ): 

/ 2 _ g-2\ M 
Q(Y,Y;qz m ) i3 = - 9 \ - -— JJ(%-j/fc X q~ 2 - y^VkQ 2 )- 



fc=i 



A. 3 Duality and the analogue of the Moore-Read state 

The p + ip model was derived from the anyonic pairing model in a limit where the pa- 
rameters a, (3 go to zero while keeping their ratio fixed, which gives the coupling constant 
G = 2/(af3~ 1 — 2) (recall (l2l]) ). One may ask wether the results obtained in this section 

4 In the context of the s-wave model this property has been exploited in |19[ 120] to study quenching 
dynamics. 
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have an analogue in the anyonic pairing model. To answer this question we first look for 
a solution of the Bethe ansatz equations flTTJ with the property that all the roots vanish, 
i.e. yj — > (j = 1, . . . , M). From eq. (TTTT) this implies 

M _2 2 2 2 

exp(4z/3(L - M - ^ + 1)) = lim TT ? 2 ~ V2 ™ = 1, M. 

The RHS of this equation is identically equal to 1 in the limit /3 — > 0, which is why there 
is such a solution in the p + ip model. One can readily check that this result remains true 
for any f3 provided 



which implies 



aioce 2 ™'", m = l,...,M, 



tf=L-£ + l = L-ff->. 



This equation coincides with the condition f )89|) for the existence of a MR state in the 
p + ip model. In the latter equation the parameter G should not be confused with the 
coupling constant of the anyonic pairing model given by Gap = sin(2/3)/ sin(a — 2/3), 
except in the limit a, (3 — > 0. The corresponding eigenstate of the anyonic pairing model 
has zero energy and it is an anyonic deformation of the MR state. Explicitly the state is 

dl 



\amr)=[J2t.) 1°) 



J=l J 

As in the p + ip model we can look for a set of roots Y, where a subset Z of P roots 
vanish, while the other set Y' with M' roots do not. The Bethe ansatz equations for the 
roots in Z yield 

M _2 2 2 2 

exp(4*/3(L - P - 2M> - ^ + 1)) = lim \ \ m " \% y m G Z 

while for the roots in Y' the Bethe ansatz equations become 

L 999 M 1 4—2 2 

±-q 2 y™4 _ -q i-rVm Vj 



exp(-2za){{ = || _ y m G F . 



Hence the roots in Y' satisfy the same Bethe ansatz equations as the original model, 
provided 

P = L-2M'- — + 1 = L- 2M' - G- 1 

zp 

which is the same condition as found in fl86l) . Hence we expect that the duality symmetry 
of the p + ip model to also be a property of the anyonic pairing model. 
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B Calculations for the p + ip model 



B.l Solution of the gap and chemical potential equations 

In this appendix we collate some results for the dimensionless parameters a and 6 which 
are solutions of equations (153(1541) . In all instances the results for the strong pairing phase 
can be deduced from those of the weak pairing phase through use of the duality relation 
( 1501) . Below, the convention for the elliptic integrals are as in [27] : 



E((f), k) = dx \J 1 — k 2 sin 2 x, 
Jo 

r<t> 

F((/),k) = dx 
Jo 



yl — k 2 sin 2 x 

With the density functions as given by (j51|52p . performing the integrals in ( I5"3"f5l| one 
finds the following results. 

ID case 

Weak coupling BCS phase: (a, b = e =F iS) 

- = VR{F{4>, k) - 2E(<f>, k)) + 2^i^±^, 
9 l + R 

R = Ve 2 + 5 2 , = 2 arctan(l/v / ^), k 2 = - ( 1 + ^) . 

2 V RJ 



where 



Moore-Read line: (a = b < 0, a = —\a\) 

1 



1 — y \a\ arctan(l/ y \a\). 



9mr 

Weak pairing phase: (a = — \a\, b = — \b\, \a\ > \b\) 



1 l + la| 



2 ' 2g 



— yfa[ E(if}, m), 



2(x-^ + ^l = J\b\F(if),m) 



where 



if; = axcsin(l/\/l + 161), m 2 = 1 - |^[. 

\ a \ 

Read-Green line: (a = —\a\ < 0,6 = 0) 



9rg 



a — v a 
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2D case 

Weak coupling BCS: (a, b = e =F i8) 



1 VTTW^2i-R + -elo, n-e + y^+W^ 



9 V -e + R 



2[x 1 = Rio, 



2 2g J °V -e + R 

where 

R= Ve 2 + 5 2 . 
Moore-Read line: (a = b < 0, a = —\a\) 

— = 1- |a|log(l + l/|a|). 
9mr 

Weak pairing phase: (a = — \a\,b = —\b\, \a\ > \b\) 

i = ^(1 + |a|)(l + \b\) - J\Ub\ + (\a\ + |6|)log f — 



x 



yi + n + v / h 7 R 

W |a6| log 1 



2 2 V V V Vl^+V^ 

Read-Green line: (a = —\a\ < 0,6 = 0) 



1 1 



1 + a — a log — — + 4/1 



fl'flG Vvl a l V '° 

B.2 Ground-state wavefunction: exact results versus mean-field 
results 

We first recall from eq. (1281) that for the p + ip model 

keK + y 



while in ID this simplifies to 

fc>0 



In the latter formulae we extended the domain of the momentum variables to include 
negative values of k x which is done to properly account for the antisymmetry of the 
wavefunction. Now we will compare the exact wavefunction with the mean-field one. 
Recalling ( 14"U|) we have 
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M 



This expression can be viewed as an average of the exact wavefunction ( 1281) which involves 
a product of pair operators characterized by the roots y m . Read and Green have shown 
1531 that in the limit where k — y the BCS wavefunction behaves as 



0(k) 



k x — iky, n < 0, (strong — coupling), 
l/(k x + iky), fi > 0, (weak — coupling). 



This behaviour can be compared with the wavefunction associated to the roots appearing 
in d2HD, 



_ k x - ik y f k x — ik y , if |k| << y/\y 

Q[ ,V) ~ k 2 -^ l/(k X + lky), if|k|»V^ 



which in real space leads to 



/oo i>oo - 
/ d 2 k^ r Q (k,y)~ mzW^y-) 
-oo J — oo I * I 

where 2 G C and i^i(z) is a Bessel function. Using the asymptotic results 
1/z (\z\ « 1) and K 1 (z) ~ e~ z ^Jix/2z (\z\ » 1) one finds 



fl(r,J/) 



z |-3/2 e -|z|V=i/ > | z | » i/yq^f, 

1/2, jzj « i/Vfy|- 



Hence at short distances all the pairs behave as in the Moore-Read wavefunction and at 
large distances they show the typical behaviour associated to localized BCS pairs with a 
correlation length 



f ^R^)' (118) 

If most of the roots y lie near the origin then the overall behaviour of the exact wavefunc- 
tion will be well described by the Moore-Read wavefunction. 

In ID the real-space wavefunction associated to a root y is given by 

/°° ■ k 
dk e ikx — = i7rsign(x) e~ |x|v/= ^ 
-oo k — y 

which decays exponentially with a correlation length given by (IllSp . 
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B.3 Wavefunction scalar product and one-point correlation func- 
tions 

The wavefunction scalar product and one-point correlation functions of the p + ip model 
can be obtained directly as a limiting case of the results in Appendix IA.2I The elements 
of the L-operator (I107P act locally, to leading order in 7, as 

L\{x) ~ l-2 P1 ^^(I-N), 
x — x~ l 

x — X 1 

Li(x) ~ ~ b , 

x — X 1 

Ll(x) ~ 7-2^^^. 

x — X -1 



Now using (11151) we have to leading order 



L 



4xp7 



z j u j 



x 2 — z 2 



0=1 3 



We define the operator (cf. ( 1251T 



^) = E^I (119) 

j=i i 

such that the states of the system are of the form 

M 

\<t>(Y)) = HC( yj )\0) (120) 
3=1 

where Y = {yj} are a solution set of the Bethe ansatz equations ( 1231) . The formula for 
the scalar products of the states (11201) follows from taking the appropriate limit of ( I112p . 
For the limit of the matrix f ll 13j) the result is 



Recalling the Bethe ansatz equations in the form (12 Op . we can substitute into the above 
to eliminate t: 

M 

F(W,Y)ij ~ AprfyjWi Yliyjwr 1 -yj x wi) 



X 



M o 



A P y ^ + 8pY 
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To simplify the calculation we introduce the diagonal matrices 



U(W,Yhj = Sij - 1 2 2 



such that 



M / M 

,...-2 



det(n(^,y» = (^j (n 
det(r(w,y)) = [p )T(w;y) 

j=1 V Vj J 



where 



T(W,Y) 



TM I _1 _1 N 



Note that 



Y(Y, Y) = C(Y,Y) 
with C(W, Y) given by f)114p . Now defining 

G{W,Y) = hmU(W,Y)F(W,Y)T(W,Y) 

leads to 

C(WY) W-^ fy ^ of w m 

W - ^ 2 ) ^ " W - 2 *) S {wl ~ wf) ^ ~ wl) 

We will also need 



X 



M M (v 2 -w 2 ) 2 

l\j} k (y]-yl)(^-^ k y 



Keeping in mind that C(x) differs from the leading term in the expansion of T 2 (x) by a 
scale factor, we now have 



n A/ n M (v 2 -w 2 ) 

rii<fc(2/i 
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Finally we can make explicit the square of the normalised wavefunction scalar product 
T(W, Y) (equivalent to the square of the fidelity) between two states as a function of the 
Bethe roots, but independent of the coupling G. The result is 



T{W, Yf 



C 2 (W,Y) det(G(W,Y)) 2 
T 2 (W, Y) det(G(W, W)) det(G(y, Y)) 

jjjj (y]-™l) 2 det(G(W,Y)f 



(y) - yl)(w] - wl) det(G(W, W)) det(G(Y, Y)) 



and we note 



M 



2 



G(Y,Y) " = %(v?-vir 

G(Y, Y)ij = {y 2_ y 2 r Hi- (121) 

In a similar manner we obtain the formula for the one-point correlation functions from 
the limit of f 1 1 1 7 j) . The result is 

I -\j \ . det(G(Y,Y) - z m W rn (Y)) 

det(G(YY)) [LZZ) 

with G(Y, Y) given by (1T2Tj) and 

B.4 Two-point correlation functions 

Finally we turn our attention to the more technically demanding problem of calculating 
two-point correlation functions. We will derive formulae for both the off-diagonal case 
(biffin) and the diagonal case (N m N n ). For the calculation of these results we extend the 
methods used in [181 E2l [69] which were developed for the s-wave model. 

First we consider the off-diagonal case. For any fixed n we may express the state 
\<j>{Y)) as 

M M 

\4>Q0) = U C ^\°) = II (p n {yp) + A£bl) |0), (124) 



p=i 



where 



c n (yp) — — l —, A 13 — — - 
n y$- Z l v$ - Z n 



Keeping in mind that (b n ) 2 = 0, b n \0) = 0, we find 

M M 

b n \<p(Y)) = J2 A ^ n )U d -(ya)\o) 

13=1 a^P 
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M M 

= J2 A nIl( c (y<*)- A M\o) 

M M M M M 

= E^n 11^)1°)- EE A n A n h i \[ ^)|0>- (125) 

p=l a^fi /3=1 l¥=a>P 

With the help of ( I125P we have that for m ^ n the unnornormalised two-point off-diagonal 
correlation functions are 

M M 

(<f>(Y)\blb n \<f>(Y)) = J2 A n^( Y )\blJ[C(y a )\0) 

0=1 a^(3 

MM M 

-EE^w y )iw n c m\°) 

13=1 a ^/3 y^a,/3 
M M 

= E^^( y )i^n^)i°) 

13=1 a^p 

MM M 

-2EE^«^)itt n c m°) ( i26 ) 

(3=1 a</3 -y^a,/3 

where in the last step we use the fact that the terms in the double sum are symmetric 
with respect to interchange of the labels a and 0. Now the two-point correlation function 
has been reduced to a sum of one-point form factors of b^ m and a double sum over the 
two-point form factors of b^b^. Recalling the definition flllOj) . for the p + ip model the 
inverse problem is easily solved as 

bl = lim y —^C(y). 

Using this expression and the wavefunction scalar product formula we obtain the matrix 
elements of both 6^ and b^ n . For the one-point case the result is 

M 2 2 M 

(<t>{Y)\bl\\C{y a M = Hm V -^(4>(Y)\C(v)l[C(y a )\0) 

a=(3 m a^P 

= X*det(72) 
where Tj^ is the M x M matrix with elements 

(7^. = Gij, j^p, 
(72) = 

*L = (yj ~z 2 J,G = G(Y, Y) is given by (HUD, and W m = W m {Y) is given by mB- (To 
keep the notation compact, we hereafter omit the Y dependence on all such matrices.) 
For the two-point form factors we have 

M 

(4>(Y)\bib{ n c( Vl )\o) 

= lim lim ^ - Z ^ v2 ~ Z H ^(Y)\C{u)C{v) f[ C(y,)\0) 

u^z m v^z n Z m Z n 

= ^ldet(T^) (127) 
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where is the M x M matrix with the elements 



( T a/3\ 
\ mnJij 


C 


3 ± a, (3, 


(rpad\ 
\ rnnJia 


= (W„) fa> 




V-mnJiB 







and 



{yl-yl){zl-zl) 



This leads to 
where 



M 



M M 



13=1 



(3=1 a<(3 



J, 



mn 



Z n(yfj — z m) 

y 2 ~ Z n 
2A a K a/3 



r /3 



2^(|/«-^)(^-<) 



We can conveniently express the determinants of T£ and in a vector notation. Defin- 
ing the M-dimensional vectors Gj , j = 1 , . . . , M and W m , m = 1 , . . . , M to have entries 

(Gj)j Gjj 

(W m ). t = (W m ) lJ 

we can write 

det(7^) = G\, Gg-i, VV m , G/3+i, Gm , 

det(T^) = Gi, Gq,_i, VV„, G a -n, ...Gg-i-, VV m , G^+i, Gm • 
Using standard properties of determinants we have 



det(T^) 



J 1/3 



Jmn 



7 (/3-2)/3 _ 



■J, 



M 



/ j u mn^^ v \ mn) "mn 

a<(3 



J 1/3 

G" mn ft fi 

i T^/T^a, <j- 

Jmn 



20 ••■■> ^ p-2 



(/3-2)/3 _^ 

— G/3-i, VV n , W m , Gm 



r 

'Jmn 



<Jmn 



Combining these we have a single determinant expression 



- mn) 



det(T^)-E^det(T°: 

Gi — z^tG 2 , Gg_ 2 ( Gg-i, Gg-i — J^ n 1 ^W n , W m , Gy3 + i, G M 



t-2/3 
On 



J, 
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This allows us to express the off-diagonal two-point correlation functions as 

(<P(Y)\blb n \<p(Y)) 



(blbn) 



(<t>(Y)\<t>(Y)) 



1 M 
det(G) Y^ J ™ det ( V " 



where 



(r)P \ - r - (^I+i yj)(yj z m) r ■ , o 1 

[ mn)i i ~ 13 (y]-yj)(y] + ,-z 2 J lU+lh 31 ' 



I mnii(^l) 



2 Z n(yfl-l Z m)(yf) Z n) 

" {yU-y}){z 2 m -zlM-zir 



Next we consider the computation of the diagonal two-point correlation functions. To 
facilitate this, we first note the commutation relation 

[N m , C(v)\ = -^bl. (128) 

With the help of the commutation relation ( 11281) . and the fact that each N m vanishes on 
the vacuum, we can write for m ^ n 

M 



{<f>{Y)\N m N n \<f>{y)) = (<f>(Y)\N m N n Y[C(y p )\0} 



13=1 

M M 

Z. 



E-Tlb-^Wm 11^)10) 

0=1 V P Zm 



MM M 

Zm \ ^ Z 



EAEA(^)i« II c <y*)\ Q ) 

MM 

a2i E 7 ,2 _ 7 2 „2 _ 7 2^ et ^") 
0=1 yP Z m V* Z n 



M M 



Z, 

=1 cl^0 



Using similar techniques as before, we can reduce the above to a sum of M determinants. 
Doing this leads to the final result 

(<j>(Y)\N m N n \<j>(Y)) 
{ m n) (<f>(Y)\<f>(Y)) 

'M-l 



V J Mp det(D p )+ J {M - 1)M det(E 

2det(G) Y J n rnn) ' J mn ^H^B 
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where 



(jJ2 i _ y|)(y 2 - z 2 ) 



i^mn) if} 







- z 2 ) 2 ' 










Gij 


(y]+i 














(y? 


- Z 2 ) 2 ' 



(E m n)ij = Gij - 2 2 — j ^ M - 1, M, 

J (y] -y M )(yj+i- z l) 

(-B'mn)i(M-l) 

(E mn ) iM 
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